A numerical relativity approach to the initial value problem in asymptotically Anti-de 
Sitter spacetime for plasma thermalization — an ADM formulation 
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This article studies a numerical relativity approach to the initial value problem in Anti-de Sit- 
ter spacetime relevant for dual non-equilibrium evolution of strongly coupled non-Abelian plasma 
undergoing Bjorken expansion. In order to use initial conditions for the metric obtained in 
[arXiv : 0906 . 4423 we introduce new, ADM formalism-based scheme for numerical integration of 
Einstein's equations with negative cosmological constant. The key novel element of this approach is 
the choice of lapse function vanishing at fixed radial position, enabling, if needed, efficient horizon 
excision. Various physical aspects of the gauge theory thermalization process in this setup have been 
outlined in our companion article arXiv: 1103.3452, In this work we focus on the gravitational side 
of the problem and present full technical details of our setup. We discuss in particular the ADM 
formalism, the explicit form of initial states, the boundary conditions for the metric on the inner 
and outer edges of the simulation domain, the relation between boundary and bulk notions of time, 
the procedure to extract the gauge theory energy-momentum tensor and non-equilibrium apparent 
horizon entropy, as well as the choice of point for freezing the lapse. Finally, we comment on various 
features of the initial profiles we consider. 
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I. INTRODUCTION 

The most fascinating theoretical challenge in the 
physics of heavy ion collisions is the understanding of 
the mechanisms behind the short - less than 1 fm/c [T] 
time after which a hydrodynamic description is necces- 
sary in order to describe experimental data. Due to the 
fact that hydrodynamics typically assumes local thermal 
equilibrium, this problem has been dubbed 'the early 
thermalization problem'. A possibility which has not 
been fully appreciated, is that viscous hydrodynamic de- 
scription may be applicable even with quite sizable pres- 
sure anisotropy |2H1] ■ Thus the quark-gluon plasma may 
still not be in a real thermalized state even though all 
the experimentally observed consequences of a (viscous) 
hydrodynamic description may still hold. 

Because of this experimental motivation, and the 
adopted usage in the heavy-ion community, we will con- 
tinue to call the transition to viscous hydrodynamics as 
'thermalization' even if the plasma is not really strictly 
thermalized in the statistical mechanics sense. The ques- 
tion, to what extent is the plasma isotropic at this 'ther- 
malization' is at the forefront of the current investiga- 
tions. 

The very low observed shear viscosity of the plasma 
produced in relativistic heavy-ion collision strongly sug- 
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gests that the plasma system is strongly coupled, thus 
making the theoretical analysis of the thermalization 
problem very difficult if not impossiblqM 

On the other hand, this observation opens up the possi- 
bihty of using methods of the AdS/CFT correspondence 
to study similar problems in gauge theories which possess 
a gravity dual. The utility of the AdS/CFT correspon- 
dence in the non-perturbative regime lies in the fact that 
then the lightest dual degrees of freedom at strong cou- 
pling are just the (super)gravity modes whose dynamics 
is given by Einstein's equations (possibly with specific 
matter fields). Moreover, from the metric sector one can 
extract the whole dynamics of the expectation values of 
the gauge theory energy-momentum tensor which is the 
key observable of interest in all hydrodynamic models. 

Thus Einstein's equations (more precisely supergravity 
equations) are an ideal arena to study the physics of ther- 
malization in a strongly coupled gauge theory. This setup 
takes the simplest form in A/" = 4 super Yang-Mills the- 
ory (SYM), when we do not excite any other expectation 
values apart from the energy-momentum tensor. Then 
the gravitational description reduces to 5-dimensional 
Einstein's equations with negative cosmological constant. 
Even more interestingly, this description is universal in 
this sector for all holographic 1+3 dimensional conformal 
field theories as argued in [6 . 

A static system in thermal equilibrium is decribed on 
the dual side by a planar black hole. It is known since 



^ In QCD we may even have a mixture of perturbative and non- 
perturbative effects. 



[7] that thermalization of generic small disturbances of 
the system is exponentially fast and is described, on the 
gravitational side, by quasi-normal modes. 

For an expanding plasma system (e.g. in the boost- 
invariant case considered here), the asymptotic equi- 
librium geometry is time dependent and looks like a 
boosted black hole and the thermalization of small dis- 
turbances is also very fast but now quasi-exponential (i.e. 
~ exp(— const r 3 )) 8J. In that paper, thermalization 
was suggested to occur as an approach to an attractor 
geometry from generic initial conditions, but clearly a 
linearized analysis is insufficicnir]- a full nonlinear treat- 
ment of Einstein's equations is needed thus neccesitating 
a numerical approach. 

Intuitively, the most important feature of gravity back- 
grounds dual to collective states of matter is the presence 
of the (event) horizon, which acts as a membrane absorb- 
ing gravitational radiation and matter outside until local 
equilibrium (in the sense of perfect fluid hydrodynamic 
description) is reached on the boundary. This mecha- 
nism of equilibration turns out to be very effective. In- 
deed numerical simulations of \TUl [TT] give short ther- 
malization times (which can be argued to correspond to 
times shorter than 1 fm/c at RHIC energies). In addi- 
tion, the result of our investigation [3] shows that viscous 
hydrodynamics applies for Trth > 0.6 — 0.7 which is con- 
sistent with RHIC assumptions (e.g. T = 500MeV and 
Tth — 0.25/m gives Tt — 0.63). These results provide 
very strong motivation for further investigations of the 
thermalization processes in the gauge-gravity dualitjrl 

What makes the thermalization process difficult are 
the many scales involved, so that the full microscopic 
description is needed. This is in stark contrast with 
the near-equilibrium dynamics, where the evolution of 
the system is governed by the equations of hydrodynam- 
ics and the only microscopic input needed are thermo- 
dynamic relations and lowest transport coefficients (up 
to first or second order in gradients). A beautiful holo- 
graphic manifestation of this fact is the fluid-gravity du- 
ality |15j , where the velocity and temperature profiles on 
the boundary specify completely the dual gravity back- 
ground and Einstein's equations can be recast in the form 
of the equations of hydrodynamics. In other words, once 
holography provides the thermodynamic and transport 
properties of the dual field theory, in order to solve the 
initial value problem in the near-equilibrium regime, one 
does not need to solve the full Einstein's equations. The 
same information is contained in the much easier equa- 
tions of hydrodynamics. 

On the other hand, in the far-from-equilibrium regime 
in a holographic field theory, one needs to specify not a 
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couple (as in hydrodynamics), but an infinite number of 
boundary functions (a couple of numbers at each con- 
stant radius slice of the initial time hypersurface in the 
bulk). Gravity then is just a clever way of recasting com- 
plicated interactions between these degrees of freedom in 
the form of classical equations of motion for the five- 
dimensional metric. 

From that perspective, the study of far-from- 
equilibrium physics leading to thermalization consists of 
two natural steps. First, one needs to specify the initial 
data for the evolution on some initial bulk time hyper- 
surface and later use fully-fledged numerical relativity 
to evolve it till the dual stress tensor can be described 
by hydrodynamics. The thermalization time is then the 
boundary time which elapsed between these two events. 

Because the thermalization process is so complicated, 
one of the most important questions is whether any reg- 
ularities emerge from underlying microscopic (here grav- 
itational) dynamics. Another way of looking at this is 
to search for some characteristic features of holographic, 
and so driven by strong coupling, thermalization, which 
might provide clues for singling out mechanisms rele- 
vant for a rapid approach to local equilibrium at RHIC. 
Yet another perspective is to try to understand various 
quantities describing thermalization at strong coupling in 
terms of some primary and secondary features of initial 
far-from-equilibrium state encoded geometrically. 

Motivated by these questions we are considering a 
boost-invariant pji- thermalization process in the holo- 
graphic conformal setting. The reason for focusing on 
the boost-invariant flow is that it is phenomenologically 
relevant for the mid-rapidity region of heavy ion colli- 
sions, but at the same time simple enough to allow for a 
thorough understanding using the gauge-gravity duality. 
In the boost-invariant case with no transverse dynam- 
ics the evolution of finite energy density system depends 
on a single coordinate - proper time r - and necessary 
starts with some far-from-equilibrium state at early time, 
which thermalizes at some transient time leaving viscous 
hydrodynamics at late time. 

The holographic studies of the boost-invariant flow 
have a long history and started with unravelling perfect 
p7] and further first [18] , second [19] and third order [20] 
hydrodynamics in the late time expansion of dual space- 
time, subsequently embedded within the framework of 
the fluid-gravity duality [15] as its special case [2T1I23) . 

The key later development underlying the present ar- 
ticle is [24], where a Taylor expansion in the Fefferman- 
Graham (FG) coordinates was used to unravel the struc- 
ture of the metric coefficients at early time and the nature 
of dual far-from-equilibrium initial state. Complemen- 
tary to the studies of boost-invariant holographic ther- 
malization reported in |10| this does not require turning 
on sources for any field theory operators and allows to ob- 
serve unforced relaxation towards local equilibrium start- 
ing from a family of far-from-equilibrium initial states. 
The crucial simplification occurring in the Fefferman- 
Graham coordinates when the initial time hypersurface 



in the bulk is chosen to be tfg = 0, is that the regularity 
of bulk geometry forces time derivatives of all 3 nontrivial 
metric coefficients (warp-factors) to vanish. This allows 
to solve explicitly 2 constraint equations and parametrize 
the most general initial data in terms of a single function 
of radial direction subject to regularity constraint and 
asymptotic AdS boundary condition. 

The coefficients of near-boundary expansion of the 
warp-factor specifying initial data turns out to be related 
to the early time expansion of the stress tensor in a one- 
to-one fashion, in accord with expectation that far-from- 
equilibrium dynamics involves many independent scales. 
By reconstructing early time power series from the near- 
boundary behavior of initial data one could not however 
see the transition to hydrodynamics, as the resulting ex- 
pressions have too small radius of convergence in t, which 
directly motivates this work. 

The task of this article is to provide a numerical frame- 
work which allows to solve Einstein's equations starting 
from initial data found in [24] and evolve them till hy- 
drodynamic regime on the boundary is reached. The 
most important physical results about the thermaliza- 
tion process in this setting were summarized in the letter 
[3] , whereas this companion article focuses on the gravity 
side of our approach and its numerical formulation. 

The principal complication in formulating a numerical 
framework is the well known diffeomorphism-invariance 
of General Relativity. The form of the equations depend 
on the type of coordinate system that one uses, and in 
addition, depend on the choice of dynamical variables. A 
good choice should be well adapted to the physical prob- 
lem under investigation and of course should be stable 
numerically. 

One edge of the computational domain can be natu- 
rally taken as the boundary of asymptotically AdS space- 
time with boundary conditions dictated by imposing the 
flatness of the metric in which the dual field theory lives. 
This condition is imposed because we do not want to 
deform the gauge theory through giving a source to its 
stress tensor or any other operator, but rather study un- 
forced relaxation from a set of far-from-equilibrium initial 
data. 

Typically one has also to impose boundary conditions 
at the other (outer) edge of the integration domain, where 
generically the curvature may be quite high. In a stan- 
dard formulation, causality of the bulk spacetime requires 
putting this edge behind the event horizon, so that ini- 
tial data specified on an initial time hypersurface encode 
boundary dynamics up to an arbitrary late time. How- 
ever, the event horizon is a global concept and cannot be 
located on a given slice of constant time foliation until the 
full spacetime is known. As a measure whether a given 
point is outside or inside the event horizon one can use 
the notion of an apparent horizon, which is intrinsic to 
a given slice. On the other hand, the existence and the 
location of an apparent horizon depends on a foliation 
and a constant time foliation induced by a given choice 
of coordinate frame might not see (at least initially) any 



apparent horizon. 

The most commonly used method of integrating Ein- 
stein's equations in asymptotically AdS spacetimes was 
to use a characteristic formulation and use ingoing 
Eddington-Finkelstein (EF) coordinates [lliniin]. This 
method has been very successful in numerical simulations 
of black hole formation in AdS spacetimes. Here, the 
outer boundary can be set to be at the location of the 
apparent horizon on the null hypersurface. Due to the 
null character no explicit boundary conditions need to be 
set there. 

For our purposes we did not adopt this formulation as 
the ingoing null light rays, do not align along tfg = hy- 
persurface and so initial conditions derived in |24| and re- 
viewed in section ITT] cannot be used as a starting point for 
the Eddington-Finkelstein approach. There seems to be 
also an unresolved technical problem in writing numer- 
ical codes in Eddington-Finkelstein coordinates starting 
from r = at the boundary stemming from the fact that 
in these coordinates the limits tef -^ and r — > oo (go- 
ing to the boundary) do not commute (see section \n\ for 
more details). 

These difficulties forced us to search for a new coor- 
dinates frame in which the hypersurface T„eui = coin- 
cides with the hypersurface tfg = and which provides 
a sensible radial cutoff in the bulk enabling numerical 
treatment. Fortunately, it turned out that successful 
code can be achieved by adopting to AdS gravity the 
ADM (Arnowitt-Deser-Misner) formalism of general rel- 
ativity |25j with fixed time hypersurfaces being spacelike. 
This scheme is also the most popular in numerical simula- 
tions in asymptotically fiat spacetimes, but so far, to our 
knowledge, has not been used in the AdS/CFT context. 
Its advantage is also that it is a very generic formula- 
tion of an initial value problem, thus it should be easy to 
generalize to other setups. 

The ADM formulation for an asymptotically AdS 
spacetime involved two difficulties which were not present 
in conventional asymptotically flat formulations. Firstly, 
the boundary conditions at the AdS boundary, dictated 
e.g. by the choice that the gauge theory metric is 
Minkowski, turned out to be surprisingly subtle and in- 
volved a careful treatment of a possible boundary diffeo- 
morphism. Secondly, and this is the chief novel feature 
of our formulation, we adopted the outer boundary con- 
ditions by freezing the evolution at the outer edge by 
making the lapse vanish there. This construction has 
several attractive features. It works perfectly well even 
when the geometry is highly curved at the edge. The ex- 
terior of the simulation domain is causally disconnected 
from the interior and thus the obtained results are com- 
pletely determined by the initial data. This last feature 
is not dependent on the location of the edge w.r.t. the 
event horizon. Thus we may perform numerically consis- 
tent simulations without any knowledge on the location 
of the event horizon. 

Let us finally note, that very recently a third type of 
numerical relativity formulation was adapted to Asymp- 



totically AdS spacetimes — the Generalized Harmonic 
formulation j26) . It would be very interesting to inves- 
tigate its relative merits with the ADM formulation (in 
particular in its more refined versions like BSSN [27J) in 
the present context. 

The plan of the article is the following. Section |ll] 
introduces the boost-invariant flow in the context of the 
gauge-gravity duality and reviews the results of [Mj . Sec- 



tion III explains how to bypass the Fefferman-Graham co- 
ordinate frame singularity by a choice of chart inspired by 
the Kruskal-Szekeres extension of Schwarzschild metric. 
Section |IV| reviews ADM formalism of general relativity, 
which in section [V| is tailored to describe the gravity dual 
to boost-invariant flow. Section|Vl]explains the subtleties 
of imposing asymptotically AdS boundary conditions and 
obtaining expectation value of dual stress tensor oper- 
ator, as well as elaborates on non-equilibrium entropy 
defined on the apparent horizon. Section |VII| describes 
various aspects of numerical side of the project. Section 



VIII elaborates on the analyzed initial conditions, com- 
pares numerical predictions for the effective temperature 
as a function of proper time with early time power series 
for three representative initial profiles, as well as explains 
subtlety in defining thermalization time for a class of ini- 
tial conditions we considered. Finally, the last section 
summarizes results and discusses open problems. 



II. BOOST-INVARIANT FLOW AND 
HOLOGRAPHY 



done on the system. It is also worth noting at this point 
that proper time - rapidity coordinates are curvilinear 
and despite the fact that the boost-invariant dynamics 
depends only on a single timelike coordinate, there is a 
hydrodynamic tail at late time in contrast to spatially 
uniform isotropization [21 [9] . 

The field theory observable, which is of interest here, 
is the expectation value of the energy-momentum tensor 
operator. This object carries direct information whether 
the system is in local equilibrium and, in holographic 
conformal field theories, undergoes decoupled dynamics 
specifying by itself dual gravity background [6j . The most 
general traceless and conserved stress tensor obeying the 
symmetries of the problem in the coordinates (IT]) takes 
the form 



T% = diag{~£,pL,PT,PT) 



(3) 



where pl and px are longitudinal and transverse pres- 
sures expressed in terms of energy density e(T) and read 



PL 



T £ 



and 



1 



Pt 



-re . 



(4) 



The precise form of the energy density as a function of 
time e(T) depends on the initial state and is governed by 
the complicated dynamics of a gauge theory or, here, by 
the dual gravity picture. The principal aim of the present 
investigation is to devise a method to obtain eir) for any 
given initial conditions. 



A. Kinematics 

Boost-invariant dynamics describes an expansion of 
plasma with an additional assumption that physics re- 
mains the same in all reference frames boosted along the 
expansion axis (i.e. in the longitudinal plane). This sym- 
metry can be made manifest by introducing proper time 
T and (spacetime-) rapidity y coordinates related to the 
usual lab frame time x*^ and position along the expansion 
axis x^ by 



= T coshy and x — t sinhy. 



(1) 



In the following x^'^ are taken to be Cartesian coordi- 
nates in the transverse plane and are denoted collectively 
as a;j_. In the absence of transverse dynamics, which is 
the simplifying assumption adopted here, the evolution of 
the system in proper time - rapidity coordinates depends 
only on proper time, since boosts along x^ direction shift 
rapidity. The background Minkowski metric in these co- 
ordinates becomes proper time-dependent 



^^boundary 



— Tj^ydx^dx^ = —dr 



'dy' 



dx\ 



(2) 



so that, as anticipated in the introduction, the system, 
by construction, is not translationally invariant in proper 
time and its evolution naturally splits into the early, tran- 
sient and late time dynamics even if no external work is 



B. Bjorken hydrodynamics and the criterium for 
thermalization 

At sufficiently late times boost-invariant plasma 
evolves according to the equations of hydrodynamics J16l . 
These are conservation equations of the stress tensor (Is]) 
under the assumption that it can be written in hydro- 
dynamic form, i.e. expressed in terms of a local tem- 
perature, velocity and gradients of velocity. Let us note 
that this assumption is not true in general. Whether one 
can write T^^, in such a form is a question of dynamics. 
Below, following [3] we will adopt an unambiguous crite- 
rion testing whether in the boost-invariant setup T^j^ can 
indeed be written in hydrodynamic form or not. 

As symmetries - boost-invariance, invariance under re- 
fiections in rapidity, as well as rotational and transla- 
tional symmetries in the transverse plane - fix the form 
of the local velocity (its only non-zero component is 
u^ ~ 1), the only non-trivial dynamical hydrodynamic 
field in the setup of interest is the local (called here 
effective for the reasons explained below) temperature 
Teff{T) defined by 



£(r) = -A^,Vr,/;(r)^ 



(5) 



The scaling of the energy density ^ with local temper- 
ature is fixed by the scale invariance of gauge theory of 



interest, whereas the prefactor counting the number of 
degrees of freedom in thermodynamic equihbrium is that 
of A/" = 4 super Yang-Mills at large N^ and strong cou- 
pling. 

We will use the equation (jsl) also in the far-from- 
equilibrium regime, but there we will treat it as a def- 
inition of an effective temperature. Physically this is the 
temperature of a thermal state of A/" = 4 SYM theory 
with the same energy density as £(t). This gives us a 
measure of the energy density which factors out the num- 
ber of degrees of freedom relevant for the specific gauge 
theory. In the rest of the text we will refer to Te// {t) as 
to an effective temperature. 

It is easy to see, that the equations of hydrodynamics 
in the boost-invariant setup reduce to a first order ordi- 
nary differential equation for the effective temperature, 
which can be solved perturbatively in the large proper 
time expansion. The result is known explicitly up to the 
third order in gradients ior Af — A super Yang-Mills and 
other (3-|-l)-dimensional holographic conformal field the- 
ories [I8H2O] reading 



A 



-21 + 27r2 + 51 log 2-24 log 2^ 



},(6) 



19447r3 {Arf 

where A is an integration constant with a dimension of 
energy governing the asymptotic scaling of the effective 
temperature with proper time [31 [TD] . Given the effective 
temperature as a function of time, A can obtained by 
fitting the hydrodynamic expression for Teff{T) to late 
time data obtained from numerical simulation. 

Introducing a dimensionless variable w being the prod- 
uct of the effective temperature and proper time 



w{t) = T^ff{T)T 



(7) 



allows to rewrite the equation of boost-invariant hydro- 
dynamics in a particularly simple form reading 



T d Fhydro{w) 

— — w — — . 

w dr w ' 



(8) 



where Fhydro{w) is a universal function of w. Hydrody- 
namic gradient expansion coincides with the laige-w ex- 
pansion of Fhydro{w). This approach to boost-invariant 
hydrodynamics is similar in spirit to previous attempts 
of Shuryak and Lublinsky of introducing all-order re- 
summed hydrodynamics [251 US] • The important differ- 
ence is that Ffiydroiw) contains all nonlinear hydrody- 
namic effects. Using the results of 18 20] we calculated 
in [3] Fhydro{w) up to third order in gradients obtaining 

Fhydrojw) ^2 _J_ l-log2 
w 3 Qnw 27'it'^w'^ 

15-27r2-451og2-h241og^2 
"^ 9727r3u;3 +■■• (9) 

Let us note that even without knowing the precise form 
of Fhydro, equation ([8]) can be used as a test whether 



T^y is in a hydrodynamic form (i.e. written purely in 
terms of a local temperature and gradients of velocity) . 
Indeed, plotting the left hand side of (Is]) as a function of 
w for various initial conditions would give a single curve if 
the hydrodynamic description would be valid or multiple 
curves in the opposite case. This analysis was performed 
in the companion article [3 . 

The equations of hydrodynamics written in the form 
([8]) allow for a simple criterion for thermalization by mea- 
suring dimensionless deviation of ([7]) from obeying ((8|. 
In particular, in [3] we adopted the following criterium 
for thermalization 



Jp3^'^ order 
hydro 



(w) 



< 0.005. 



(10) 



Note that the condition ( 10 1 is based on demanding that 



the effective temperature obeys equations of hydrody- 
namics, rather than on isotropy of the pressures. Indeed, 
as the results of [5] and earlier studies in [TO] show, the 
pressure anisotropy can be quite sizable, nevertheless the 
evolution of the system being governed by hydrodynam- 
ics. In section [VIII| we discuss the sensitivity of thermal- 
ization times obtained from ( 10 ) to the number on the 



right hand side of ( 10 1 



Once gradient terms become negligible, the entropy is 
no longer changing in time. This can be seen by evalu- 
ating the entropy per unit rapidity and transverse area, 
being a product of the thermodynamic entropy 



s{r) = ^N^n^T^ffirf 



(11) 



and the volume element scaling linearly with proper time 
(pi). In the following, as in [3], in order to measure the 
final entropy, we will get rid of N^ factor as in ^ and will 
be using instead dimensionless quantity, being entropy 
per unit rapidity and transverse area measured in units 
of the initial effective temperature T^ 



eff 



dS/dydx\ 



(12) 



Applying this definition in the r — >■ 00 limit of the hy- 
drodynamic regime gives us an expression for the final 
(dimensionless) entropy expressed in terms of A 



A^ • (7i;^) 



(13) 



In the latter part of the article we will be using a partic- 
ular generalization of entropy to non-equilibrium regime 



which however reduces to (131 in the regime of 



•^n—eq^ 

applicability of perfect fluid hydrodynamics 

C. Early time dynamics from holography 

In the gauge-gravity duality the symmetries of the 
boundary dynamics are also the symmetries of the dual 



background, which here is a solution Einstein's equations 
with negative cosniological constant 



1 d{d-l) _ 

Rab — ^Rgab 7^ ~ ^ 



(14) 



with d being the dimension of boundary spacetime taken 
here to be 1 + 3. In the following we set L - the radius of 
AdS vacuum solution - to 1. The most general (4 + 1)- 
dimensional dual metric sharing the symmetries of the 
boost-invariant flow takes the form 

gabdx^'dx^ = -{- A^df + t^B^dy^ + 



"'°bulk 



1 



C^dx^ + —D^du^ + 2Edtdu], 

-^ All J 



Au 



(15) 



where the warp-factors are functions of bulk proper time 
t and radial coordinate u with u — denoting the bound- 
ary [30 . Note that t does not need to coincide with r even 
at the boundary. The use of u instead of z = y'u soft- 
ens singularities of Einstein's equations at z = 0, which 
is convenient when solving them numerically. There is 
a redundancy in the metric (151 coming from diffeomor- 



phisms in t and u, which can be used to fix two out of 
three warp-factors within A, D and E. Various choices 
define various coordinate frames covering different parts 
of the underlying manifold in various ways. 

The choice D — and E = ~7r7= leads to coordinates 
of ingoing Eddington-Finkelstein type used in numeri- 
cal simulations of [H (TUl [H] (typically the Eddington- 
Finkelstein coordinates are written using r — l/^/u). 
One reason why in this work we are not using the 
Eddington-Finkelstein coordinates is that there seems to 
be a subtlety there in starting at r = at the boundary 
(being also t = in the bulk). This can be seen by looking 
at the vacuum AdS metric in the Eddington-Finkelstein 
coordinates 



1 



ds^ = dt^ 

u 



1 



t 



dy^ 



1 2 1 

—dx, —j-dtdu, 

u u-^l ^ 

(16) 

for which [dy^ term) near-boundary limit [u — > 0) does 

not commute with small time limit t — > 0. The authors of 

[lOj have not encountered this problem as their numerical 

simulations were always starting at tim = Tmi > 0. 

Setting D = 1 and E = in the metric Ansatz ( 15 ) 

leads to the Fefferman- Graham chart, which has been 

used extensively to understand the asymptotic properties 

of bulk spacetime relevant for obtaining the stress tensor 

using the procedure of holographic renormalization [31] . 

In particular, the boundary condition setting the metric 

in which dual field theory lives to be flat is just 



ApG = 1, BpG = 1 and Cpc = 1 



(17) 



and in that case one can also identify bulk t and bound- 
ary (physical) proper time t. In the FefFerman-Graham 
coordinates the near-boundary expansion of warp-factors 
turns out to occur in integer powers of u. The first sub- 
leading term vanishes automatically and the first non- 
trivial term in the expansion of Apc upon holographic 



renormalization [3T] comes out proportional to the energy 
density of dual Af = 4 super Yang-Mills plasma £(t) 

TT^ 9 n^ A , 1 „, Q 
^--^-iVf^"-4iVf(/+3^>+-- 

BFG-l-f,{s + re'y - ^y + irs">^ + . . . 



(18) 



Nc in the formulas above denotes the number of colors, 
which although formally infinite, cancels out with N^ 
contribution in the energy density yielding in the end a fi- 
nite result. All the remaining terms in the near-boundary 
expansion of warp-factors turn out to be entirely speci- 
fied in terms of the energy density and its derivatives. 
It is worth noting that one way of making sure that 
dual spacetime written in other coordinate frames has a 
fiat boundary (i.e. imposing such boundary condition), 
as well as obtaining energy density of dual field theory 
configuration, is to perform a near-boundary coordinate 
change to the Fefferman-Graham chart and then directly 
use ( 17 1 and ( 18 1 to identify the relevant terms. This will 



become important later in the article. 

The starting point of the analysis of the early time dy- 
namics in [2j was the expansion (18). Because of the 



appearance of inverse powers of proper time in terms 
containing odd derivatives of energy density, demand- 
ing finiteness of near-boundary warp-factors in the limit 
T — )■ gave a highly nontrivial constraint on the early 
time Taylor series of energy density, so that only even 
powers of proper time are allowed 



eir)\ 



= ^N^TT^ 



T 



(») 



eff 



^ rpll 

2^eff 



{oy 



(19) 



One consequence of this is that first proper time deriva- 
tives of warp-factors at r = vanish, so that the two 
gravitational constraint equations (the uu and tu compo- 
nents of Einstein's equations) provide solvable relations 
between AiT-G (r = 0,u) = Ao{u), Bpg{t — Q,u) = Bo{u) 
and Cfg(j = 0, m) = Co(m) and their radial derivatives 
(since r derivatives vanish at t = 0). One of these rela- 
tion can be explicitly solved giving 



Bo{u)=A^{u), 
whereas the other leads to a differential equations 
Kin) , C'^{u) 



Aq{u) Co(u) 



0. 



(20) 



(21) 



This equation, via taking suitable combination of Aq i 
Co can be solved exactly (i.e. in terms of a single func- 
tion specifying initial condition, but with no integration 
involved) , as shown in ^24, , but for the purposes of this 



article, it will be sufBcient to solve (21) numerically for 
Aq regarding Cq as an initial condition. This way of 
viewing ([2T|) has the advantage that by passing to an- 



other coordinate frame in which initial time hypersur- 
face coincides with FefFerman- Graham one, Cq will not 
transform (up to redefinitions of what is meant by u). 
An interesting feature of the equation (21 ) is that for all 
allowed Co{u) there exists a point uq such that Ao(u) 
has a single zero there [Mj . This point signals the break- 
down of the Fefferman-Graham chart and makes it hard, 
if not impossible, to evolve Einstein's equations in this 
coordinate frame. A typical, analytic example of initial 
conditions derived in [5T is 



Aq = C0S7U and Cq = cosh 7 m. 



(22) 



where 7 



^V3TT^{T^ff) sets the initial energy density 
7r/27 with the Fefferman- 



'eff) 



{FG) 



and u runs from to u, 
Graham coordinate frame breaking down at the latter 
point. Table H] in the appendix summarizes 29 different 
Co{u) and corresponding values of Uq considered in 
this article. 

Although the near-boundary expansion of Cq's stays in 
an one-to-one correspondence with the early time Taylor 
series for the energy density ( 19 1 



Co{u) 



l-'iT^%r-' 



n^(T^^ 



'rjill 



{0)u^ 



(23) 
and given Cq in an analytic form (at least close to the 
boundary) one can work out the first couple of dozen 
terms in the expansion of s{t), the radius of convergence 
of the resulting series is insufficient t o obs erve the transi- 
tion to hydrodynamics (see section VIII for an explicit 



comparison of the effective temperature given by the 
early time power series with the full numerical result). 
Not surprisingly, in order to understand thermalization 
in this model, a full numerical solution of the initial value 
problem on the gravity side is needed. 

A naive attempt at numerically solving Einstein's 
equations in the Fefferman-Graham frame is bound to fail 
for two reasons. Firstly, the generic vanishing of Ao{u) 
leads to a coordinate singularity. Secondly, even if this 
was overcome, some sensible initial conditions exhibit a 
curvature singularity in the bullrl(at u — >■ 00). 

The first step towards solving numerically Einstein's 
equations is thus to find a coordinate frame, which allows 
to use initial data found in [21], bypasses the breakdown 
of Fefferman-Graham coordinates and makes it possible 
to introduce bulk cutoff for numerical simulation in order 
to avoid (a possible) curvature singularity. 



III. NEW COORDINATES 

The way we deal here with the singularity of Fefferman- 
Graham frame at t = is analogous to what happens in 
the case of Schwarzschild black hole 



ds^ = -(1 )dt^ + (1 ; 

r r 



Hr'^+r'^dn% (24) 



when going from Schwarzschild to Kruskal-Szekeres co- 
ordinates. The latter ones are defined by 



V^^(^_l)i/V/4^^sinh-47, 
^ 2M ' AM ' 

C/ = (^^-l)i/V/4^^cosh — 
^ 2M ' AM 

and lead to manifestly regular metric at r = 2M 



ds^ = ^^^e~'^/2^^(-dy2 + d[/2) + rHn 

r 



(25) 
(26) 

(27) 



with r defined implicitly by the relation V^ — U^ = (1 — 
_r_-|gr/2M^ What is important for our purposes is that 
the hypersurface i = coincides with the one given by 
T^ = 0, as follows from ( [25| . In that way we can take t — 
metric functions at dftg2 (analogues of dy^ and x"^ warp 
factors) and use them directly as dO^a (in our case dy"^ 
and x^) warp-factors setting initial data for the evolution 
in a chart without spurious coordinate singularities (here 
the analogue of Kruskal-Szekeres coordinates). To make 
the analogy even sharper, on the V — hypersurface 
one can use r instead of U coordinate at the price that 
Schwarzschild metric will become IZ-dependent. 

Let's try to apply the same logic to the boost-invariant 
metric in the Fefferman-Graham coordinates 



1 t^ 

ds^ = AFG{t,u)df H BFG{t,u)dy'^ 

u u 

H — CFG{t,u)dx'^, + --^du^ 
u Aw' 



(28) 



in the neighborhood oi t — hypersurface. The coordi- 
nate singularity arises from Apcit = 0, m) = Bpcit = 
0, u) vanishing at some radial position, whose precise 
value depends on the choice of initial condition. By doing 
a local coordinate transformation in the neighborhood of 
t = (i.e. perturbatively in t) one can redefine i in a 
fashion similar to ( 25 ) and give an arbitrary form to the 
dt^ metric coefficient at this particular instance of tim^ 
This statement pesists to an arbitrary order in small t 
expansion and can be understood as adopting a different 
gauge within the choices allowed by the metric Ansatz 



* Such initial conditions may be physical if there is an event hori- 
zon cloaking the curvature singularity. This wrill turn out to be 
the case for the initial conditions considered in the present work. 



An explicit example of such a redefinition is 






i+o{P) 



( 15 1 - one in which one fixes A and E at the cost of leav- 



ing D dynamical. This is advantageous for us, as we can 
use now ADM formalism-based scheme for numerically 
integrating Einstein's equations. Moreover, the freedom 
of choice of A allows to introduce a very natural bulk 
cutoff as anticipated in the introduction and elaborated 
on in the next section. Last, but not least, t = hyper- 
surfaces in both coordinate frames are by definition the 
same and one can use u as a radial coordinate on both 
of them. This allows us to start with initial conditions 
derived earlier in the Fefferman- Graham coordinates and 
continue the early time power series solution of [24] be- 
yond its radius of convergence in order to explore the 
transition to hydrodynamics, which was one motivation 
for the present work. 



IV. REVIEW OF ADM FORMULATION OF 
GENERAL RELATIVITY 

A particular formulation of Einstein's equations which 
is very convenient for studying evolution from generic 
initial data is the ADM formulation [531 123 • In this 
work we did not adopt any of its more refined versions 
like BSNN [27], as in 1-hl (and even 2-hl) ordinary ADM 
should suffice. 

The key idea behind the ADM formalism, making it at 
the same time a natural point of departure in implement- 
ing Einstein's equations numerically, is to assume that 
there exists a global foliation of spacetime into spacelike 
hypersurfaces of constant time and recasting Einstein's 
equations in terms of equations intrinsic to a constant 
time slice (constraint equations) and extrinsic ones (en- 
coding the time-evolution) [32]. As it will turn out, due 
to the choice of coordinates we will be using foliations of 
patches of the bulk spacetime, rather than global folia- 
tions. 

Denoting by A a scalar function slicing the bulk onto 
constant time hypersurfaces, the induced metric describ- 
ing the intrinsic geometry of a leaf takes the form 



lab = gab +nanb, 



(29) 



where n° is a future directed unit normal vector obtained 
from the gradient of A. The induced metric (29) acts 



also as the projector operator onto foliation leaves. The 
spacetime embedding of the constant time hypersurface 
is described by the extrinsic curvature 



Kr, 



1 



~.*^7i^ ab-) 



(30) 



where £„ denotes the Lie derivative along the normal 
direction n°' . A coordinate basis temporal vector, can be 
constructed from the unit normal n° as 






= t" = 5n° 



/3" 



(31) 



with a and /?" being respectively the lapse and the shift 
vector {rf'Pa = 0). The role of the lapse is to measure 



the rate of time flow between subsequent slices of the 
foliation, whereas the shift vector describes how the hy- 
persurfaces slide onto each other in transverse directions. 
With the use of the projector (29), Einstein's equations 



Rab — ^Rgab — SttGat Tab 



(32) 



can be decomposed into constraints and evolution equa- 
tions. The bulk energy-momentum tensor (not to be con- 
fused with the expectation value of the boundary stress 
tensor operator!) is decomposed into density, current and 
a transverse tensor taking respectively the form 

p = Tabn''n\ j, =. -Tabu'^-il S^d = Tabl\l\. (33) 

In the ADM formulation, equations governing internal 
spacelike geometry of the hypersurface are recast in 
the form of constraint equations reflecting the time and 
space decomposition of spacetime. The Hamiltonian con- 
straint, following from Gauss equation, reads 



R + K^- KabK 



ah 



IQttGnP, 



(34) 



whereas the momentum constraint derived from Codazzi 
equation is 



D,K''„ 



DaK = SnGNJa. 



(35) 



Da here is the covariant derivative compatible with the 
spatial metric jab ^^nd TZ is the Ricci scalar associated 
with it. Evolution equations for the induced metric come 
from projecting its Lie derivative onto constant time slice 



Ctlab 



-2aKab+DaPb + Dt/3a- 



(36) 



The evolution equations for the external curvature can 
be obtained in a similar fashion from the Ricci equation 



CtKab = -DaDba + aiUab - 2KacK'b + KabK) 

+ P'D.Kab + K.bDaP'' + K,aDbP' 

p — s 



8TrGNa[Sab 



d-l 



lab\ 



(37) 



Here TZab is the Ricci tensor of 7afc, S — S°'a ^^^d d denotes 
dimensionality of boundary (taken everywhere in the pa- 
per, apart from this section where it is kept general, to 
be 1 -|- 3) . One can notice that the only place where the 
bulk stress tensor contributes is the last term and this is 
where vacuum AdS cosmological constant will manifest 
itself. Gomparing ( 14 ) with ( 32 ) one can see that in the 



absence of matter fields the bulk energy-momentum ten- 
sor Tab is related to the radius of vacuum AdS solution 
L, {d + 1) -dimensional Newton's constant Gjv and bulk 
metric gab through 



Tab — 



d{d-l) 



(38) 



By introducing transverse coordinates of the foliation leaf 
y' one can recast the metric into the standard ADM form 



ds' 



a^dy + j,j{dy' -f- I3'd\){dy^ -f pUX). (39) 



By doing so one can confine (d+l)-diniensional indices to 
d spatial slices as in the projected quantities normal com- 
ponents are zero (we shall denote those by i and j Latin 
letters, as opposed to a, b and c symbols running through 
d + 1 indices). Moreover, now Lie derivative along n" is 
simply a time derivative dx (since A parametrizes curves 
normal to slices). 



V. ADM FORMULATION FOR 
BOOST-INVARIANT PLASMA 

Before we apply the ADM equations in the context of 
the time evolution of the boost-invariant geometry in an 
asymptotically AdS geometry, we have to specify certain 
additional ingredients. 

Firstly, due to the fact that the metric blows up at the 
AdS boundary, one has to redefine the ADM variables i.e. 
the spatial metric coefficients and the extrinsic curvature 
by taking out predefined factors which will ensure that 
the asymptotic behaviour is fulfilled while keeping all the 
new redefined variables finite throughout the integration 
domain. 

Secondly, we have to specify the initial conditions 
which satisfy constraint equations. In the special case of 
a boost-invariant geometry with initial conditions set at 
t = 0, this requires some care as the initial hypersurface is 
not strictly spacelike but has signature (0, -I-, -I-, -I-). This 
feature simplifies the determination of possible consistent 
initial conditions, but also requires special treatment of 
the first step of the temporal evolution. Of course, for 
all t > 0, the constant t hypersurfaces are spacelike and 
conventional ADM formulation applies. 

Thirdly, we have to impose boundary conditions for all 
dynamical variables at the AdS boundary. These bound- 
ary conditions are conceptually clear, as they correspond 
to enforcing a Minkowski metric on the boundary (in or- 
der to ensure that dual A/" = 4 SYM theory is defined 
on an ordinary Minkowski space). However, it turns out 
that for a generic choice of lapse functions, the boundary 
metric may be related to Minkowski by a boundary diffeo- 
morphism. Taking this into account leads to more com- 
plicated boundary conditions than conventional Dirichlet 
boundary conditions. 

Fourthly, we have to impose boundary conditions at 
the outer edge of the integration domain in the bulk. 
These boundary conditions are much more subtle and are 
not fixed by the very principles of AdS/CFT (as are the 
previous boundary conditions) but rather depend on the 
specific features of the dynamical problem at hand. The 
main requirements are that these boundary conditions 
should not interfere with the physics of interest and also 
should lead to stable numerical evolution. 

Finally, we have to specify the final ingredients of the 
ADM formalism - the lapse and shift functions which 
encode how the hypersurfaces of fixed (simulation-)time 
fit together to form the 5D geometry. For our purposes 
we set the shift vector to zero, however the lapse will 



remain nontrivial. A part of the lapse function will be 
used in defining the outer edge boundary conditions but 
the remaining part will have to be specified. 

In the following we will discuss in turn all the above 
mentioned points. 



A. Rescaled ADM variables and equations of 
motion 



In order to ensure that the functions entering the ADM 
equations are finite everywhere, even at the AdS bound- 
ary, we have factored out appropriate factors of u (the 
AdS bulk coordinate). Our (4 + 1) -dimensional metric 
takes the form 



2_ -a^{u)a^{t,u)dt^ t^a^{u)h'^{t,u)dy'^ 



ds^ = 



+ 



&'{t,u)dx\ d^{t,u)du^ 



4u2 



(40) 



with empty AdS being represented by all the coefficient 
functions equal to unity. Note that in general the time 
coordinate t will not be equal to the gauge theory proper 
time r at the boundary. We will derive an explicit rela- 
tion between the two coordinates in subsection lV CI With 
the above definition, the ADM spatial metric jij is given 
by 



7ij — diag 



t^a'^{u)b'^{t,u) c^it,u) c^{t,u) d'^{t,u) 
u ^ u ' u ^ Au^ 



(41) 

The nontrivial components of the extrinsic curvature are 
also rescaled 



K, 



diag 



ta{u)L{t,u) M{t,u) M{t,u) P{t,u) 



Au\/u 



Finally, the lapse function is parametrized by 

a{u)a(t,u) 



(42) 



(43) 



The reason of factoring out the time independent func- 
tion a(u) will be clear when we discuss the outer bound- 
ary conditions below. For simplicity we will always set 
a(0) ^ 1. 

In the presence of a cosmological constant A = —6 and 
vanishing shift the vacuum ADM equations become 






(44) 



dtK, 



2a{u)a(t, u) 



a{u)a{t,u) a{u)a{t,u) 

1= I 4 -1= 'JijV'i '^j: 



{R,j-2K'^K,,+KK,j] 



For the metric coefficients we obtain: 

db{t,u) _ -b{t,uf + a{t,u)L{t,u) 
dt ^ th{t,u) '' 

dc{t,u) a{u)a{t,u)M{t,u) 
dt " c{t, u) ' 

dd{t,u) a{u)a{t,u)P{t,u) 



(45) 



dt 



d{t,u) 



The evolution equations for the extrinsic curvature func- 
tions (3tL, dtM, dtP are quite lengthy and can be found 
in Appendix |B] They were generated by Mathematica 
and directly transfered to the computer code. 

The hamiltonian and momentum constraints take the 
form 






(46) 



with only the u component of the momentum constraint 
C4 being a-priori nontrivial. 



choices made, the final initial conditions for the ADM 
evolution are 



bo{u) = do{u)= Loiu) = 1 Moiu) = Po{u) - 



co(u) 



profile 
^0 



{u) 



(49) 



parametrized by the single function c^^°-'^ '^{u). Com- 
paring the resulting initial geometry with the Fefferman- 
Graham initial condition we see that we can identify 
^ro/i e^^-j ^j-|-]-^ -j-j^g Fefferman-Graham initial condition 
Cfg{t = 0,u) discussed in subsection II C Conse- 
quently, we have at our disposal a power series solution 
for e(r), which may be used to check the results of the 
numerical evolution for some initial range of t. 

In some cases, when running the simulations, we no- 
ticed that quite narrow structures emerge at late times 
close to the outer edge of the simulation domain in the 
bulk. This causes the numerical evolution to eventually 
break down for a given size of the spatial grid. In these 
cases we found it useful to redefine the initial coordinates 

by 



B. Initial conditions at t = and the first step of 
the numerical evolution 

Typically in the ADM formalism, the initial conditions 
are obtained by solving the constraint equations ( 46 1 . In 



the case of the boost-invariant geometry, however, the 
initial hypersurface i = is not spacelike as it contains 
the light cone in the longitudinal plane. Its signature is 
(0, -I-, -|-, -|-). This requires an ab-initio analysis of ini- 
tial conditions, which become in fact unconstrained, and 
a special treatment of the evolution equations (i.e. the 
right hand side of dtb, dtc etc.) at t = 0. 

To this end we will expand the lapse and all coefficient 
functions up to linear order in t: 



b ^ bo{u) + bi{u)t - 

d — do{u) + di {u)t - 

M = Mq{u) + Mi{u)t- 

a = ao{u) + ai{u)t - 



c — co(u) + ci{u)t+ . . 
L ^ Lq{u) + Li{u)t + . 
P ^ Poiu) + Pi{u)t + 



(47) 



After inserting these expansions into both the ADM dy- 
namical equations (44) and the constraints (46), we ob- 



tain the initial conditions (feoiCo etc.) as well as evolu- 
tion equations at t — 0. We found by an explicit calcu- 
lation that both co(u) and do{u) are unconstrained and 
free, whereas fog (u) turns out to be proportional to ag (u) . 
Without loss of generality we will set the constant of pro- 
portionality to unity. All this taken into account gives 

bo{u) = ao(u) Loiu) = ao{u) Mq{u) = Pq{u) = 0. 

(48) 
Note that we are free to perform a spatial diffeomorphism 
(redefine u). In this way we may freely set do{u) = 1, 
leaving the initial condition to be completely specified by 
a single function Co(u). In the following we will restrict 
our lapse functions to satisfy ao(u) = 1. With these 



1-Ci 



(50) 



with appropriately choosen constant C. This redefinition 
stretches the grid at the outer edge allowing in some cases 
for longer evolution. The initial conditions now take the 
form 



do{u 



boiu) ^ Loiu) ^ 1 Moiu) = Poiu) = 0. (51) 



Of course all the physics extracted from running the sim- 
ulation from the initial conditions (49) and (51 ) with the 



same function c^^°^^ '^(m) is completely equivalent. 

The terms linear in i in (47) give the right hand sides 



of the evolution equations at i = 0, which we use for the 
first time-step of the numerical integration. We obtain 
in this way 



9tfo(0,u) = 0, 9tc(0,u)=0, 

dtdiO,u)^0, atL(0,0)=0, 

-2ai2uicduciaa) + duiaacduc)) 
dtM(0,u) = 

aadudic? + duC^ 



(52) 



d2 

2aac? 2aac^ 



d3 
dtPiQ,u) ^ 'iudliaa) + 



u ud^ 

Adudiaa + duiaa)) 



2aa(l — d^) Auaad^c Auaaducdud 
u c cd 

As a final note let us clarify why in the present ADM 
formalism we have completely unconstrained initial data, 
while in the Fefferman-Graham case we had the differen- 



tial constraint (21). The Fefferman-Graham coordinates 



are a special case of the general metric ansatz (40) al- 
beit with the constraint that (i(i, u) = 1 (imposing this 
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condition, however, transforms the lapse into a dynami- 
cal variable). Requiring that the condition d{t,u) = 1 is 
preserved under evolution requires that P{t, u) = Q and 
in particular dtP{0, u) = 0. It is this last equation which 
reduces exactly to the the Fefferman-Graham initial data 
constraint (21 ). 



C. Boundary conditions at the AdS boundary 

at It = 

The boundary conditions at u = are choosen as to 
ensure that the gauge theory metric is just the (3 + 1)- 
dimensional flat Minkowski metric. The easiest way to 
derive these conditions is to start from the leading asymp- 
totics in Fefferman-Graham coordinates (i.e. basically 
the empty AdS^ geometry) 



ds^ 



-dr 



FG 



' FG 



dy^ + dx\ + dz^ 



(53) 



and consider the most general change of coordinates (in 
an expansion around u = 0) to our ADM metric ansatz 



( 40 1 . Hence we write 

TFG = k{t)+h{t)u + f2{t)u^ + . 

z = go{t)u^ + gi{t)ui + g2{t)ui 



(54) 



The physical gauge theory proper time is just r = /o(i). 
Inserting this into ( 53 ) we see at once from the dv? com- 



ponent that the boundary condition for d is just 

d(i,0) = l. (55) 

Consequently P{t,Q) = 0. However the boundary values 
of b and c can be t-dependent. We get that go{t) = 
l/c(i, 0) (from dx\) and combining this with 



fo 
50 



= a{0)b{t,0)t 



(56) 



(from dy^) we obtain a very important relation between 
the coordinate time t in our ADM simulation and the 
physical proper time (recall that we have a(0) = 1) 



fo{t) 



b{t,0)t 



(57) 



The condition that ensures having the Minkowski met- 
ric will be the compatibility of the above equation for 
/o(t) with last remaining condition coming from the dt^ 
component of the metric 



/o = 



a(t,0) 
c(t,0)' 



(58) 



This condition leads to an expression for the boundary 
value of L 



L(t,0) =b{t,0) + t 



c2(i,0) 



M{t,0). 



(59) 



For actual numerical implementation it is convenient to 
impose the above boundary condition in a differential 
form 



dtL = dt [b + t-M 



(60) 



At this stage we are missing only the evolution equation 
for M{t, 0). There are various ways of implementing this 
equation. Perhaps the simplest way is to expand the 
evolution equation for dtM around u = and take into 
account the conditions derived above i.e. d(t,0) = 1, 
P{t, 0) = and the expression for L{t, 0). 

An alternative way is to compute the first subleading 
terms in the coordinate transformation ( |54[ ). In any case 
this computation has to be done in order to derive the 
formula for the energy density £(t). Requiring that there 
is no dudt term, one can show that 



hit) 



9ait)ga{t) 
Voit) 



(61) 



gi (t) can be computed by looking at the first subleading 
term in the dt^ component of the metric and can be ex- 
pressed in terms of a, dua, 9„a, as well as /o(i), <?o(0 ^^'^ 
their derivatives. Finally, we may use these expressions 
to compute 9„d (at u — 0) from the first subleading term 
in the du^ component. After taking into account all the 
above definitions and relations, the result may be brought 
to the form (all expressions are evaluated at u = 0) 

dud = -2dua h -^-^ 5-. (62) 

This equation allows us to solve for dtM which is the last 
remaining equation at the AdS boundary. 

To summarize, the boundary conditions at u = are 
formulated as the following evolution equations for the 
boundary values of the ADM variables 

-52 + aL 



tb 



aM 



dtb{t,0) = 

dtc{t,0) 
dtd{t,0) = 0, 

dtL{t,Q)^dt{b + t^M), 

dtM{t, 0) = -2c2(a„(aa) + adud) + 
dtP{t,0)^0. 



3aM 



/f2 



2 2c2 



(63) 



As a final point let us comment on the regularity of the 
ADM evolution equations in the bulk as we approach 
M — > 0. The evolution equations for the metric coef- 
ficients are explicitly regular. The right hand side of 
the evolution equations for the extrinsic curvature coef- 
ficients have, however, the following structure 

dt{L, M, P}^ - (d^ - 1) ■ finite + regular. (64) 
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Since our boundary condition for d is d{t,0) = 1, the 
above term does not lead to any numerical problems. 
This is partly because we use spectral methods which 
preserve very well the smoothness of functions close to 
the boundary. 



D. Boundary conditions at the outer edge of the 
simulation domain 



As explained in section II, there is a subset of ini- 
tial conditions for which the curvature goes to infinity in 
the center of AdS. These initial conditions still may be 
perfectly physical if the singularity would be surrounded 
by an event horizon. However a-priori we do not know 
where the event horizon is located. We can locate it only 
after running a simulation. Moreover due to the null 
character of the initial hypersurface at t = 0, the condi- 
tions for an apparent (dynamical) horizon cannot be met 
there. 

In order to perform the simulation we have to cut off 
our numerical grid at some finite value oi u = ug and im- 
pose boundary conditions which would not contaminate 
the true physical evolution. The more or less standard 
choices in numerical relativity cannot be applied here. 
Putting any boundary conditions inside the event hori- 
zon cannot be used in our case as the location of the 
event horizon is unknown when we start the simulation. 
Moreover, the spectral discretization that we use makes 
the setup very sensitive to the boundary conditions since 
Chebyshev differentiation is very much nonlocal in con- 
trast to finite difference discretization. The second stan- 
dard choice, namely imposing outgoing radiative condi- 
tions at the outer edge of the simulation domain is also 
not an option, as the geometry at the outer edge may be 
highly curved. 

In order to bypass the above difficulties, we have de- 
cided to use the freedom of the choice of spatial foliations 
in the ADM formalism by requiring that all hypersur- 
faces (i.e. for any fixed t) pass through the same single 
spacetime point u = Uq on the initial hypersurface t = 0. 
This can be done technically by freezing the evolution 
at u = Wo which amounts to forcing the lapse to vanish 
there. We achieve this by setting the t-independent part 
of the lapse in ( 40 ) to be 



A) 



a(u) 



cos 



2 Uq 



(65) 



All ADM evolution equations are regular at u = uq and 
indeed do not lead to any instabilities. Moreover, the 
region of spacetime outside our numerical grid is causally 
disconnected from the simulation domain. Namely, the 
asymptotic boundary of our simulation domain will be a 
null geodesic running from the spacetime point u — uq 
and i = in the direction of the boundary. 

However we must note that the optimal choice of uq is 
crucial c.f. Fig. IT] If uq lies between the event horizon 
and the boundary, the null geodesic from uq will reach 
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FIG. 1. Schematic view of constant time foliations depending 
on the locus at which the lapse vanishes (no). Situation A) 
shows a simulation filling in a triangle covering only a finite 
interval of boundary time. By pushing the position of bulk 
cutoff inwards (until one passes the event horizon), one can 
recover the dynamics on larger and larger regions of bound- 
ary. In the situation B) the lapse vanishes at the position of 
the event horizon on the initial time hypersurface, which en- 
sures maximally long simulation time (theoretically infinite) 
and is optimal for studying the transition to hydrodynamics 
in dual gauge theory. In the case C) the simulation pene- 
trates into the black brane/black hole interior revealing the 
apparent horizon at the cost of breaking down when constant 
time slices start approaching the curvature singularity. This 
type of behavior can be seen explicitly in Fig. [2] showing the 
results of a sample simulation. 



the boundary at some time t — t^. Then clearly the 
ADM simulation will never be able to go past t = t^, and 
the simulation will break down there. The optimal choice 
would be to set uq to be exactly at the position of the 
event horizon. Then, the simulation could in principle be 
run indefinitely and the simulation domain would fill the 
whole exterior of the event horizon. This choice is the 
one we are aiming at in order to extend the simulation 
to late times and to observe the transition to hydrody- 
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namics. Also in this case, the event horizon ensures that 
the simulation would never encounter the curvature sin- 
gularity. If, on the other hand, u ^ uq would lie beyond 
the event horizon, the simulation would break down at 
some finite time due to the curvature singularity. This 
choice of uq is in fact also useful, as it allows us to locate 
the apparent (dynamical) horizon and e.g. measure the 
entropy of the plasma system. This choice is also crucial 
in determining the optimal value of uq as explained in 
subsection IVIBI below. 

As far as we know, the outer boundary condition ob- 
tained by freezing the lapse has not been used in the 
literature so far. 



E. The choice of lapse functions 

Even after fixing a{u) we still have to determine the 
final ingredient in the ADM formulation - the remaining 
part of the lapse function a{t,u). If we were to adopt 
the simplest choice a{t,u) = 1, then we would encounter 
coordinate singularities (e.g. for the profile cq — cosh7M, 
the d{t,u) warp- factor develops a zero). This is in fact 
a very well known behaviour of the ADM formalism in 
asymptotically flat spacetime p3 35 , which can be reme- 
died by standard choices of dynamical lapse functions, 
like Y^det 7^^ or 1 -I- log det jij . However these choices 
do not seem to work well in the case at hand. We have 
chosen a couple of lapse functions by trial and error. The 
rationale was e.g. to make a proportional to d to avoid 
the vanishing of d and make it inversely proportional to 
b in order to avoid a too quick rise of b. 

We always normalized our lapse functions to be equal 
to 1 at t = 0. Therefore we set 



a(t, u) 



lapse{t,u) 
lapse{0, u) 



(66) 



The choices which seem to work best, i.e. gave long 
enough bulk evolution to see thermalization in the 
boundary theory, depending on the initial profile were 



lapsei — 



dc^ , bd 


lapses - 


d 


b ' ^"P"^' 1 + ^62' 
and lapses = d. 


b 



(67) 



Another issue which influenced the choice of lapse func- 
tions is that for some of them one could ensure that the 
simulation time t coincided with the proper time on the 
boundary. This happens when 6 = c at u = 0. In order 
to preserve this equality under evolution, it is enough to 
require that0a(t,O) = c(t, 0) (or fo(t,0)). So from the 
above lapses, lapsei and lapse2 ensure the equality t = t 
at the boundary. 



VI. OBSERVABLES 

With all the above ingredients in place one can run 
the numerical evolution, which is described in more de- 
tail in section VII. Once the geometry is known, we will 
be interested in extracting the energy-momentum tensor 
from the near-boundary behaviour of the metric. An- 
other quantity of interest will be the determination of 
the appearance and the precise location of an apparent 
horizon and the extraction of a nonequlibrium entropy of 
the plasma system. 



A. Energy density and transverse and longitudinal 
pressure 

Obtaining the components of the energy-momentum 
tensor from the results of the numerical simulation in the 
ADM variables turns out to be surprisingly subtle. The 
main complication comes from the fact that the radial 
du^ component of the metric is a nontrivial dynamical 
time-dependent field and in addition our ADM time co- 
ordinate may differ from the physical Minkowski proper 
time. 

The way to perform the computation is to perform 



the change of variables ( 54 1 from the Fefferman-Graham 
form 



ds' = 



-{I - ^-^e{T)u')dTla + r|G(l + ^Pl(t)«^) V 



(1 + ^-^pT{T)u')dxl + ^du^ 



(68) 



and compare the result with the near boundary expan- 
sion of the ADM metric (40). This requires the com- 
putation of the terms /2(i7and 32 (i) in (54 1 which is 
straightforward but requires the use of Mathematica. We 
also repeatedly use the boundary equations of motion 
(63). Then we may extract £{t), plIt) and pt{t) from 



expressions for d^d{t,0), d'^b{t,0) and d^c{t,0) respec- 
tively. Note that we extract directly p^ and px, even 
though they could be obtained from £(r) using (|4|, for 
two reasons - we avoid taking a temporal derivative of 
the numerical data for e(r) and as a byproduct can check 
whether the numerical evolution preserves the conserva- 
tion of the energy-momentum tensor. 

After carrying out the calculations mentioned above 
we arrive at the following formula for e(r) 



sir) 



'27r2 



d'a- 



2duadua dla 



duadud 



duOtdud 



\d^d^-\dld- 



dudtP 
Aa 



(69) 



This follows directly from equation J58l 



and similar formulas for pl{t) and pt{t). This for- 
mula still has a drawback, as it includes dependence on 
dudtP{t, 0). However we can evaluate this expression by 
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taking the u derivative of the equation for dtP{t, u) and 
expanding this near the boundary. In this way we eh- 
mated all time derivatives from the right hand side of 
e(r). The final expression for £{t) is 

sir) =^c^ I dta + 9„a9„d 



2ducdud 9 



dud^ 



duP iMduP 



Ub 



4c2 



The expression for pT is given by 



(70) 



Pt{t) = g hdlc - Ic^d^d' + \c^dld + \d^dM^ 



1 



c^MduP 



(71) 



The longitudinal pressure pl{t) can be found from trace- 
lessness e(T) = Pl{t) + 2pt{t). 

The above expressions are quite complex and in order 
to test them for possible errors we have compared eir) 
extracted using the above expressions from the numerical 
simulation with the power series solution for e(T) and 
found complete agreement (see Fig. [t] later in the text). 
A further test was to extract e(r) from two simulations 
which used different lapse functions. Again we found 
agreement. 



B. Apparent horizon entropy density 

Apart from the energy density and longitudinal and 
transverse pressures, a very important physical property 
of the evolving plasma system is its entropy density per 
transverse area and unit (spacetime) rapidity. Accord- 
ing to the AdS/CFT correspondence, gauge theory en- 
tropy can be reconstructed from the Bekenstein-Hawking 
entropy of a horizon. This is well established and un- 
ambiguous in the static case, however in the far-from- 
equilibrium time-dependent setup the situation is much 
more subtle. In fact, it is not even clear whether an exact 
local notion of entropy density makes sense on the QFT 
side. However, phenomenological notions of local entropy 
density are widely used in dissipative hydrodynamics. 

In the present work, we adopt a natural geometrical 
prescription for local entropy density which reduces to 
the one used in dissipative hydrodynamics in its regime 
of validity. 

On the AdS side, the definition of entropy density in- 
volves two distinct steps. In the first step, one calculates 
the area element of a (certain kind of a) horizon. The 
Qth Qj.(jgj. requirement that one has to impose is that the 
horizon area only increases i.e. the horizon satisfies an 
area theorem. Moreover, causality has to be preserved 



i.e. the horizon area cannot increase in anticipation of 
some event. This condition rules out the use of event 
horizons [21 [3S] . Currently, the most natural choice is 
the use of apparent horizons described in detail below. 

In the second step, one has to link the area element of 
the horizon to a specific point on the boundary in order 
to associate the local entropy density (area) to a definite 
spacetime point in gauge theory. We follow the proposal 
of [37] of shooting a null geodesic from the point at the 
boundary and taking the area element from the point of 
intersection of the null geodesic with the apparent hori- 
zon. 

Below we discuss these two steps in more detail. 

Apparent horizons are quasilocaQ notions of black hole 
boundaries, always confined to the causal interior of a 
black hole [3S1[5H]- From the point of view of this paper, 
they are of interest for two reasons. In the first place, 
their existence is useful, because an intrinsic ("local") 
notion of crossing the boundary of black hole can be uti- 
lized to adjust the radial cutoff for integrating Einstein's 
equations in numerical simulations. Secondly, similarly 
as for event horizons, their area only increases thus sat- 
isfying an area theorem. These properties make them a 
good ingredient for a causal generalization of entropy to 
non-equilibrium situations in holographic systems. All 
this is in stark contrast with the event horizon with its 
teleological nature and acausal evolution 39J. 

The notion of an apparent horizon, and so entropy de- 
fined on it, is tied to a particular foliation of spacetime, 
which leads to ambiguities. Indeed, in order to find an 
apparent horizon, one first defines the foliation of space- 
time into constant time slices. Then on each slice one 
locates, if it exists, such a codimension-2 hypersurface 
that wavefronts of light rays emitted (into the future) 
in an outward direction (i.e. here towards the bound- 
ary) stay constant in area, whereas wavefronts emitted 
in the inward direction shrink. The union of all these 
codimension-2 hypersurfaces, being itself a codimention- 
1 hypersurface, is what we call here, with a slight abuse 
of terminology, an apparent horizon (in fact this is the 
definition of a future outer trapping horizon). Closer 
considerations reveal that the area of constitutive slices 
of apparent horizon is never decreasing, leading to area 
theorem [39] . 

In highly symmetric spacetimes, such as gravity duals 
of boost-invariant flows, one typically searches for appar- 
ent horizons respecting the physical symmetries. Other- 
wise, field theory entropy extracted from the apparent 
horizon would not obey the symmetries of the state un- 
der consideration [40] . This assumption vastly simplifies 
searches of an apparent horizon in the gravity dual to 
the boost-invariant flow, as null directions of ingoing and 
outgoing wavefronts, denoted in the literature by null 



^ Their definition requires the existence of fully trapped surfaces in 
the neighborhood of an apparent horizon, hence the term quasilo- 
cal. 
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vectors 1°" and n", are fixed by symmetry, up to normal- 
ization factor, after imposing 



I J" = 0, nau" ^ and lau" 



-1. 



(72) 



The latter condition specifies cross-normalization of the 
vectors with the minus sign on the right hand side en- 
suring that Z° is future pointing if n" is and vice versa. 
In particular, for the metric Ansatz (40 1 vectors solving 
(ItI read 



^,^/,(t,.)|- "("y'") , 0,0,0,-%^), (73) 



1 J a{u)a{t,u) 
h(t,u) I V^^/u 



0,0,0, 



2V2u _ 

V^d{t, u) 

2\/2u 



,(74) 



where h(t, u) is a positive scalar function playing no role 
in the discussion here. 

The condition for an apparent horizon on a given time 
slice, which is taken here to be i = const, is that the rate 
of change of the area of a lightfront (transverse area of 
congruence of null geodesies) in 1°- direction is 0. This 
quantity is measured by the expansion 9i given by 



9z = (5''^ + ;V + zV)v,/, 



(75) 



with On defined in a completely analogueous way. 

In practice, for a given constant time hypersurface in 
coordinates fixed by the metric Ansatz (40) with some 
definite choice of lapse, one searches for w's such that di 
vanishes and On is negative. Fig. [2] shows the apparent 
horizon for a particular choice of lapse and sufficiently 
large radial cutoff. In all the simulations we did there 
was no apparent horizon on the initial time slice. How- 
ever, it was always possible to choose large enough radial 
cutoff to see it develop later during the simulation. Our 
coordinate frame makes it appear, starting from some 
time, at two radial positions on each constant time slice, 
forming a U-shaped structure surrounding the curvature 
singularity, as shown on Fig. [2] 

The apparent horizon's entropy density is defined here 
as Bekenstein-Hawking entropy (density) 



SAH 



—a{u)b{t, u)c{t, u) 



AGn u 



(76) 



Left as it is, it is still not a local entropy from the dual 
field theory point of view. We still need to associate 
points on the horizon with points at the boundary by 
providing a so-called bulk-boundary map. Such an asso- 
ciation is not free of ambiguities, however here it seems 
that the only sensible way to do so is to associate horizon 
points with points on the boundary lying on the same in- 
going radial null geodesic [20, 36j, much in the spirit of 
Vaidya spacetime and its generalizations [H]. Moreover, 
due to the large number of symmetries, the direction of 
the null geodesic is essentially fixed and unambiguous. 
Note that such mapping is consistent with the one intro- 
duced in the context of fluid-gravity duality [37] . 
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FIG. 2. The apparent horizon (black curve) and a radial 
null geodesic (red curve) sent from the boundary (left edge 
of the plot) at r = into the bulk for a sample profile (no. 
23 from table n|. This curve coincides with a curve of fixed 
Eddington-Finkelstein proper time tef = 0. Background col- 
ors correspond to curvature with blue denoting small cur- 
vature and red large curvature. The red region inside the 
apparent horizon denote the neighborhood of the curvature 
singularity. One can see (and check the numerical factor ex- 
plicitly) that curvature at the right edge of the plot remains 
the same during the evolution, in agreement with expectation 
that vanishing lapse freezes the time flow there. The curva- 
ture at the left edge of the plot also remains constant due to 
imposed AdS asymptotics. 



For all the considered initial states, we found a non- 
zero apparent horizon entropy at the boundary time 
T = 0, in the sense explained above. Although in the 
non-equilibrium regime there is no clear microscopic pic- 
ture of apparent horizon entropy, its initial value com- 
pared to the final one is a useful measure of how far from 
equilibrium a given initial state was. In particular, using 
that Gat = ^^ |32] one can reexpress (76) in terms of 



dimensionless entropy density Sn-eq measured in units of 



eq 



the initial effective temperature T^ A 



Sn—c 



SAH 



as introduced in section Hi Bl 

VII. THE NUMERICAL SIMULATION 
A. Determination of uq 



(77) 



As explained in section [VD[ in order to be able to ex- 
tend the simulation to sufficiently large values of proper 
time and to observe a clean and unambigous transition 
to hydrodynamics, it is crucial to optimally tune the 
position of the cutoff in the bulk which we denote by 
uq. The optimal value of mq would be the position of 
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FIG. 3. Determination of the cutoff point Uq ensuring 
long evolution time necessary to see thermalization in the 
boundary theory for a sample profile (no. 23 from table |l|. 
In order to obtain Uq , which is expected to approximate 
well the position of the event horizon at initial time hypersur- 
face, we anticipate that for large enough time the position of 
the event horizon coincides with the position of the apparent 
horizon we consider and shoot backwards in time an outgoing 
null geodesic (almost) tangent to late time apparent horizon 
(plotted as thick red curve). Other outgoing null geodesies 
are plotted as arrowed curves. The apparent horizon is plot- 
ted as thick black curve. Background colors, as in Fig. [2] 
correspond to curvature, from blue (small curvature) to red 
(large curvature). 



tion uq will be a very good estimate of the initial position 
of the event horizon and will allow for a long period of 
evolution. 

The initial exploratory simulation was also important 
for us as we used it to extract the initial entropy cor- 
respond ing to the given initial conditions as outlined in 
section (IVIBl). 



B. The numerical implementation and tests 

In the present work we adopted a Chebyshev dis- 
cretization of the spatial grid. This allows us to use 
quite modest grid sizes {N = 201, 257, 513, 1025 depend- 
ing on the initial profile) and also allows for very accu- 
rate evaluation of spatial derivatives at the AdS bound- 
ary which are necessary for extracting the gauge theory 
energy-momentum tensor. The downside is that Cheby- 
shev differentiation is very much nonlocal so any prob- 
lems at the edges will affect quicky the whole integra- 
tion domain. We implemented the Chebyshev differenti- 
ation using Discrete Cosine (Fourier) transforms and the 
f ftw3 library. Time stepping was done using an adap- 
tive 8*''/9*'' order Runge-Kutta method from the GNU 
scientific librarj|j 

In order to test the numerical simulations in all cases 
we monitored the preservation of the ADM constraints 
in the form 



the event horizon on the initial hypersurface. Then, 
apart from purely numerical complications, simulations 
could be extended indefinitely and the simulation domain 
would completely cover the whole exterior of the event 
horizon, safely separating the evolution from the curva- 
ture singularity. 

The principal stumbling block is the lack of knowledge 
of the position of the event horizon at t = 0. Moreover, 
due to the lightlike nature of a part of the initial hyper- 
surface, we cannot locate an apparent horizon there using 
the vanishing of expansion scalars. 

However, the practical determination of the optimal uq 
is in fact quite simple. First we run our simulation with 
a relatively large value of uq for which the simulation 
would break down due to the appearance of a curvature 
singularity (as shown in Fig. IC). An initial guess for 
this first exploratory value of ug is typically provided by 
taking uq to be 10 — 20% larger than the position of the 
FG coordinate singularity listed in table |l] 

The representative outcome for this exploratory sim- 
ulation is shown in Fig. [3J We identify the apparent 
horizon using the vanishing of the expansion scalar 9i. 
To get a first estimate of the position of the event hori- 
zon at i = 0, we take the outgoing radial null geodesic 
from the neighborhood of the late time outer edge of the 
apparent horizon and evolve it backwards in time until it 
reaches the initial time hypersurface. The resulting posi- 
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(78) 



and similarly for the momentum constraints. As an addi- 
tional check, for some profiles we compared the energies 
£(t) extracted from simulations done with two different 
choices of ADM lapse for the same initial conditions. We 
found a very good agreement. Another completely inde- 
pendent cross-check was to compare the e(T) extracted 
from the numerical simulation with s{t) computed ana- 
lytically as a power series in r. Within the radius of con- 
vergence of the power series, we found excellent agree- 
ment (some example comparisons will be presented in 
the following section). Finally, for a fixed profile (no. 23 
in Appendix A), we also checked whether the dynami- 
cal horizon identified for simulations with two different 
choices of lapse functions really coincide. This was done 
by i) comparing the areas extracted from the intersection 
of the dynamical horizon with null geodesies propagat- 
ing from the boundary for a range of r, ii) comparing 
the curvature invariants Rabcd.R"''""^ at these intersection 
points. Again we found perfect agreement. 
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FIG. 4. Plots of initial warp-factors as functions of square root of radial position. Thicker curves denote parts of initial data 
outside event horizon on initial time slices. Numbers in the legend match profile numbers collected in table IT] (color online). 



VIII. THE PROFILES 

In the companion article '3] , we have presented a thor- 
ough analysis of the physical characteristics of thermal- 
ization of a boost-invariant plasma system obtained using 
the numerical setup described in the present paper. 

In this section, we would like to discuss in more detail 
certain aspects which are interesting in view of the results 
obtained in p]. Firstly, we discuss various geometrical 
characteristics of the initial conditions (of which we give 
a complete list in Appendix A) in relation to the event 
and apparent horizons, then we describe the qualitative 
behaviour of the proper time evolution of the effective 
temperature giving evidence that a transient rise of the 
tempearture in the non-equilibrium regime observed for 
initial conditions with high entropy is a real effect and 
not a numerical artefact. Finally, we discuss stability 



issues of our criterion for thermalization (10). 



A. Geometric characteristics of the initial states 

We considered 29 different initial conditions specified 
by Co(w), as shown in table |l] in Appendix A. We focused 
on profiles, which had no singularity for finite valueaj 
of w, but we allowed for a possible blow up as u — >■ oo. 
Within the considered class, the numerical evolution from 
initial profiles gave rise to quite rich dynamics, exhibit- 



ing a wide variety of behaviors. Nevertheless, in [3] we 
observed surprising regularities - entropy extracted from 
the apparent horizon at r = turns out to characterize 
the key features of the transition to hydrodynamics. It 
is therefore interesting to understand this in more detail. 
As the initial time hypersurface coincides with the 
one in the Fefferman-Graham coordinates, we can use 



Fefferman-Graham constraint equation (21) and locate 



.(FG) 



radial position Uq for which Ao{u) vanishes (the co- 
ordinate singularity of Fefferman-Graham chart, as orig- 
inally pointed out in "M] and outlined in section In]) . 



Ao(4^'=))=0. 



(79) 



Although we have no clear geometric interpretation of 
this particular point, it turns out that there is a surpris- 
ing correlation between the initial non-equilibrium en- 
tropy, as obtained from an apparent horizon using the 
prescription elucidated in section VI B and the radial 



position of Fefferman-Graham singularity at i = (see 
Fig. 11). 

The position of the event horizon on the initial time 
slice i = 0, is to a very good degree approximated by 
the value of uo determined as in figure Js] We found that 
it is also correlated with the position of the Fefferman- 
Graham singularity (see figure l6| , and hence with the 
initial apparent horizon entropy. This result was very 
helpful in running the actual simulations, as it provided 
a good first guess for the optimal position of the radial 
cutoff. 



It would be interesting to relax this assumption. 
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FIG. 5. Dimensionless initial non-equilibrium entropy given 



by ( 77 1 as a function of the position of Fefferman-Graham sin- 



gularity showing clear correlation. Both quantities are plotted 
in the units of the effective temperature at t = 0. Color code 
matches figures B] and m (color online) . 



As outlined in [3], it seems that, at least for the class 
of profiles considered here, the apparent horizon entropy 
at r = provides a key infrared characteristic of the 
initial conditions. One possible explanation, supported 
by figure |4j is that the event horizon on the initial time 
slice effectively removes from the dynamics the part of 
spacetime for which there is a significant difference be- 
tween various initial data. As a result, the space of ini- 
tial data on the initial time slice truncated at the event 
horizon might be actually quite simple and crudely char- 
acterized by a couple of parameters of which the initial 
non-equilibrium entropy together with the initial effec- 
tive temperature may be the most important ones. Note 
also that the initial non-equilibrium entropy is measured 
on tef — Eddington-Finkelstein hypersurface differing 
from t = hypersurfaces where our initial conditions are 
imposed. The close correlation between Uq (on t = 0) 
and the initial entropy (defined on tef = 0) is, there- 
fore, quite surprising. It would be interesting to under- 
stand the initial value problem in the ingoing Eddington- 
Finkelstein coordinates and check whether the mecha- 
nism behind the simplicity of initial states conjectured 
here is indeed correct. 



B. Qualitative features of plasma expansion 
— cooling and reheating 
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FIG. 6. Approximate position of the event horizon on the 
initial time slice as a function of the position of Fefferman- 
Graham singularity showing clear correlation. Both quanti- 
ties are plotted in the units of the effective temperature at 
T — 0. Color code matches figures |4] and [s] (color online). 



conditions characterized by largest initial entropy among 
the states we considered, led to a growth of the effective 
temperature in the initial stage of evolution followed only 
later by a cooling phase, as required by Bjorken hydro- 
dynamics. The peculiar feature of the latter states was 
that the effective temperature at thermalization was for 
them sometimes higher than the initial one. This does 
not mean that thermalization (understood here always as 
a transition to hydrodynamics) occurred instantaneously 
but rather that it occurred after a sizable non-equilibrium 
evolution including a reheating phase. Finally, we also 
encountered a profile, with an intermediate value of ini- 
tial non-equilibrium entropy compared to others we con- 
sidered, exhibiting initially a plateau in the effective tem- 
perature as a function of time, which only later decayed. 
For the three types of the initial states discussed above 
we compared the energy density obtained from numeri- 
cal simulation with truncated early time power series ob- 
tained directly from the initial profile, as introduced in 

We observed perfect 



and reviewed in section II C 



agreement within the convergence radia of the power se- 
ries solutions. The agreement constitutes another non- 
trivial test of our numerics. This also shows that the 
reheating behavior is a real effect and is not due to some 
pathologies in the numerics. 

In all these cases we also see that the radius of conver- 
gence of the power series is smaller than the time of the 
transition to hydrodynamics. 



The investigations of the transition to hydrodynamics 
in the companion paper [3j revealed that for some initial 
conditions, the plasma does not cool but initially under- 
goes a period of 'reheating'. Fig. [7] shows the time evo- 
lution of the effective temperature for three sample pro- 
files. Depending on the value of initial non-equilibrium 
entropy, we found three types of behavior. Profiles with 
small initial entropy, as compared to the initial effective 
temperature, led to the effective temperature (and so the 
energy density) decaying quite rapidly with time. Initial 



C. Stability of the criterion of thermalization 
— subtleties for small entropy initial data 

Finally, we would like to point out a subtlety in deter- 
mining the thermalization time, especially relevant for 
profiles with low initial non-equilibrium entropy. In [5] 
we defined thermalization time as the time after which 



^■^w (with w defined as in ([7|)) obtained numerically 
from gravity deviates from third order hydrodynamics re- 
sult ^ by less than 0.5% (see the criterion (lOl). As the 
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FIG. 8. Thermalization level (red line) and criterion expres- 
sion (blue line) plotted as a function of time for the profile 
no. 29 in table ll with s^'l^^ = 0.4761. 
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FIG. 7. Comparison of effective temperatures as functions 
of time obtained from numerics (solid gray and dashed blue 
curves) and given by early time power series expression (dot- 
ted red curves) for three sample profiles. Situations A), B) 
and C) correspond respectively to initial profiles 7, 23 and 
29 from table III representing initial states with small, inter- 
mediate and large initial entropies. Solid gray curves denote 
far-from-equilibrium part from the evolution, whereas dashed 
blue curves correspond to hydrodynamic regime with ther- 
malization time set by the criterium (10 1. Each plot shows 



perfect agreement between numerical results and early time 
power series solutions within the radia of convergence of the 
latter. Note also that the radius of convergence of early time 
power series is much too small to observe the transition to 
hydrodynamics. 



magnitude of acceptable deviation from hydrodynamics 
is a free parameter (at least within a reasonable range), 
thermalization is never a clear cut event. Moreover, the 
form of Ffiydro is known only up to the third order in 
gradients and so it might happen that the thermaliza- 
tion criterium ( |10[ ) is actually sensitive to the inclusion 
of terms higher order in gradients (fourth order hydro- 
dynamics and higher), whose precise form has not been 



FIG. 9. Thermalization level (red line) and criterion expres- 
sion (blue line) plotted as a function of time for the profile 



no. 3 in table ll with s 







J^i 



- 0.0200. 



determined so far. 

As an example of this subtlety, let us present two fig- 
ures with the time dependence of the left hand side of 
our thermalization criterium. Figures |8] and |9] contain 
plots of the expression (10 1, which compares the val- 



T-j^w and 



ues of (numerical) exact function F{w) 
the one obtained from third order viscous hydrodynamics 
Phydrcf'^^ i''^) ■ Thermalization time is determined numer- 
ically from the last intersection point between the plot of 



of (10) and the threshold line at 0.005. 



In Fig. |8] the thermalization time is determined in a 
very robust way as within a reasonable range of variation 
of the threshold, there is just a single point of intersection 
with the threshold line. Moreover, the curve ( 10 ) is very 



steep so modifying the threshold value by a factor of 2 
or 4 (0.01, 0.02) would lead to only a small change in the 
determined thermalization time. 

Fig. [9] shows an analogous plot for an initial condition 
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with one of the smallest initial non-equilibrium entropy. 
In contrast to Fig. [8] we have three intersection points 
which are quite unstable w.r.t. changing the threshold 
(0.01, 0.02). Indeed, for these higher threshold values 
there would be just a single intersection point located 
closer to the previous example. The slowly decaying tail 
of ( 10 1 suggests that it may be of a hydrodynamic nature 
albeit of a higher order type (4*'' order and higher). We 
expect that including the currently unknown higher order 
hydrodynamic terms into Fhydro{w) would bring Fig 
to look more like Fig. [8] Indeed, knowing in princip. 
the all-order Fhydro{w), the function 



dr 

Ffiydroiw) 



1 



(80) 



should be exponentially small on the right hand side of 
the above plots. 

We adopted the criterion with a fixed threshold in or- 
der to avoid ambiguities related to subjective judgement 
which point of intersection to choose. However, we have 
to take into account that the small initial entropy data 
are much more subtle in this respect. 



IX. SUMMARY AND OPEN DIRECTIONS 

This paper introduces an ADM formalism-based 
scheme for numerical integration of Einstein's equations 
with negative cosmological constant in the setup describ- 
ing gravity dual to a boost-invariant strongly coupled 
plasma system. The main motivation for this came from 
the earlier work by some of us [21], which used the 
Fefferman-Graham coordinates to constrain the form of 
initial time (r = 0) metric components specifying in this 
way gravity representation of dual far-from-equilibrium 
initial states (see section In]). In [24] we also studied time 
evolution in a power series form, albeit with a too small 
radius of convergence to observe directly a transition to 
hydrodynamics. Our present numerical approach over- 
comes this shortcoming. We are interested, as in [24!, m 
unforced (i.e. with all sources in the gauge theory turned 
off) relaxation of these states in order to clearly separate 
the thermalization process from the creation of the ini- 
tial non-equilibrium states (as opposed to earlier works 

mm)- 

Numerical treatment of the initial value problem re- 
quires specifying it on a compact domain with boundary 
conditions not altering the evolution, at least outside of 
the event horizon. Unfortunately, the Fefferman-Graham 
coordinates used in ^24J do not provide a sensible no- 
tion of bulk cutoff, in particular, based on the example 
of AdS-Schwarzschild black brane they are not expected 
to cover interiors of dynamical black branes of interest. 
Because of this, we introduced a new coordinate frame, 
which, much in the spirit of the relation between Kruskal- 
Szekeres and Schwarzschild coordinates in the case of 
Schwarzschild black hole, coincides with the Fefferman- 
Graham chart only at the initial time hypersurface, but 



extends beyond the Fefferman-Graham coordinate singu- 
larity. Moreover, this new coordinate system allows for 
a convenient bulk cutoff and, if needed, for recovering 
black brane interior or horizon excision (see section III 
and VD). 

The new coordinate frame we introduced allows to 
adopt ADM formalism-based scheme for numerical evo- 
lution (see sections IV andlvl. The key technical element 
of the present work, which to our knowledge has not been 
explored in the literature before, is a new treatment of the 
outer boundary condition in a potentially highly curved 
spacetime through forcing the lapse to vanish at a fixed 
radial position. This stops the flow of time at a given 
point (actually on a codimension-2 hypersurface) leading 
to a very convenient bulk cutoff for numerical integration, 
as we are free to specify its position. 

In order to ensure a long enough simulation time to 
see the transition to a hydrodynamic regime in the dual 
gauge theory, such a cutoff should be very close to the 
position of the event horizon on the initial time slice. 
Although the latter information is not available from the 
start due to the teleological nature of the event horizon, 
we nevertheless managed to find a simple way of deter- 
mining it. To do so, we ran a trial simulation with large 
enough bulk cutoff to eventually see an apparent hori- 
zon forming and then traced back outgoing radial null 
geodesic tangent to the late time apparent horizon. As 
in the late time regime we expect our apparent horizon to 
coincide with the event horizon, the trajectory of outgo- 
ing null geodesic of interest should, to a very good degree, 
coincide with the location of the event horizon. The in- 
tersection point of this trajectory with the initial time 
hypersurface provides the locus where the lapse needs to 
vanish in order to ensure a long simulation time. 

Another nontrivial issue which we encountered, is con- 
nected with imposing asymptotically AdS boundary con- 
ditional^ as well as obtaining the expectation value of the 
gauge theory energy-momentum tensor carrying crucial 
information about the thermalization process (see sub- 
sections VC and VI A[ ). The chief difficulty came from 
the fact that due to the dynamical time-dependent ra- 
dial metric coefficient, the points on constant radius hy- 
persurfaces approach the asymptotic boundary at var- 
ious rates. This complicates both obtaining the form 
of the boundary metric, as well as invalidates the "stan- 
dard" formula for dual stress tensor expressed in terms of 
extrinsic curvature of constant radius hypersurface. To 
circumnavigate this issue we performed near-boundary 
coordinate transformation to the Fefferman-Graham co- 
ordinates, in which it is particularly simple to obtain ex- 
pressions for the boundary metric and expectation value 
of the gauge theory energy-momentum tensor. 



^^ This was an important point, as we wanted to ensure that the 
plasma system we are studying evolves in ordinary Minkowski 
spacetime and not in some time-dependent background geome- 
try. 
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One interesting finding of our studies is that, because 
of fixing bulk gauge freedom by specifying the form of a 
time-like warp-factor, the boundary and bulk notions of 
time differed by a boundary diffeomorphism. This also al- 
tered the form of the Minkowski metric on the boundary. 
We expect that similar issues might arise when trying to 
impose normalizable boundary conditions when study- 
ing e.g. spectrum of quasinormal modes in spacetimes 
expressed in coordinate frames with the gauge fixed by 
specifying the form of warp-factors associated with the 
field theory directions. 

The accuracy of our numerical simulation, based on 
spectral discretization in radial direction and finite dif- 
ference time stepping using high order Runge-Kutta 
method, has been monitored in a couple of ways (see 
section VII). In the first place, we used two appropri- 



ately normalized constraint equations, which, if small at 
all points of the grid, provide a nontrivial check of the 
numerics, as we used an unconstrained evolution scheme. 
The other way of making sure that our results are cor- 
rect, was to run the numerics for a given initial profile 
with different choices of lapse (both the position at which 
the lapse vanishes and its functional dependence on the 
dynamical metric components). As various choices of 
lapse cover various parts of underlying manifold foliating 
it in various ways by constant time hypersurfaces, ob- 
taining for each choice the same (up to numerical error) 
effective temperature as well as apparent horizon entropy 
as functions of boundary time gave highly nontrivial in- 
dications of the accuracy of numerical simulation. The 
third way of making sure the numerics worked fine was to 
compare the effective temperature as a function of time 
obtained numerically with its analytic form in the early 
time power series (see section VIIIp ). Again, within the 
radius of convergence of the early time power series we 
observed perfect agreement. 

A very surprising finding of this work, reported first in 
the companion article [3 , is that characteristics of ther- 
malization for all the profiles we considered are correlated 
with the initial value of non-equilibrium entropy defined 
on the apparent horizon (see subsection |VI B for details 
on what we mean here by local non-equilibrium entropy) . 
A possible explanation of this might be that the event 
horizon appears on the initial time hypersurface as soon 
as the initial time warp-factors start to differ significantly 
from their vacuum value. It would be very interesting 
to verify this finding by considering initial conditions in 
the Eddington-Finkelstein frame, as in this case, due to 
constant time foliation being aligned along radial ingo- 
ing null geodesies, the apparent horizon would possibly 
appear from the start (see however (2). It would be 
thus interesting to understand the correlation between 
the position of the apparent horizon (if it exists) on the 
initial time hypersurface with the level of complication 
of a warp-factor setting initial conditions. It would also 
be interesting to perform explicit numerical simulations 
and check how the position of the event horizon is cor- 
related with the position of the apparent horizon on the 



initial time hypersurfaces in the Eddington-Finkelstein 
coordinates. 

Another interesting point is related to the choice of 
initial profiles we considered here. We focus only on the 
profiles which are either regular for all radii or blow up as 
the radial coordinate is taken to infinity. In principle, it is 
possible that there are initial profiles diverging for a finite 
radial position, but having the event horizon shielding 
the singularity. It would be interesting to investigate 
this issue in more detail, as such studies might uncover 
a new class of initial states. 

The profiles we considered exhibit a wide variety of 
behaviors in the bulk. Although, as expected on general 
grounds, the evolution of the one-point function of the 
stress tensor do not depend on the details of a given pro- 
file hidden behind the event horizon, nonlocal observables 
such as entanglement entropy, Wilson loops or higher 
point correlation functions generally do [43j. It would 
be thus interesting to investigate this issue in our setup, 
as has been done for the Vaidya spacetime in [3H H5] . 
Our setup might be advantageous compared to setting 
initial conditions in Eddington-Finkelstein type of coor- 
dinates when it comes to obtaining equal time two-point 
functions of heavy scalar operators, as these are approx- 
imated by spacelike geodesies, which might closely align 
along our initial constant time slices. 

Yet another interesting direction to follow would be to 
perform linear stability analysis of the early time boost- 
invariant plasma with respect to perturbations breaking 
boost-invariance (i.e. depending on rapidity), as such 
instabilities have been found in the weak coupling regime 
within the colour glass condensate approach j46j . 

One direction, which has not been explored at all so far, 
is to understand holographic thermalization in the grav- 
ity backgrounds dual to non-conformal strongly coupled 
gauge theories. It would be of obvious interest to study in 
such gauge theories boost-invariant flow, as the simplest 
phenomenologically interesting model of plasma dynam- 
ics with the use of holographic methods introduced here 
or in [TU]. 

Finally, it would be fascinating to introduce transverse 
dynamics in the holographic model of Bjorken flow to 
get a handle on thermalization phase of the elliptic flow 
I47j , as well as understand how nontrivial structures in 
the transverse plane at the initial time affect the initial 
conditions for hydrodynamic evolution |48j. One moti- 
vation for doing this, stemming from the results of this 
work reported in [3] (and also earlier in [TU]), is that 
in the one-dimensional boost-invariant flow non-Abelian 
plasma can be very anisotropic, yet well-described by hy- 
drodynamics. Allowing the plasma system to expand in 
traverse directions should lower, at least intuitively, the 
degree of anisotropy at thermalization. Investigating this 
issue in more detail is important, as this touches upon the 
question whether the QCD quark-gluon-plasma produced 
in heavy-ion collisions indeed needs to be approximately 
isotropic before the transition to hydrodynamics. 
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Appendix A: A summary of the initial conditions and simulation results 



TABLE I: Below we collect information about the initial profiles we considered: Cq{u) is the initial condition as 



discussed in section II with 7 = ^^''^'^{'^eff) ! "0 is a position of corresponding coordinate singularity in the 

Fefferman- Graham chart; Mq is the approximate position of the event horizon on the initial time slice; w^^'^' is 
thcrmalization time measured in terms of an effective temperature at thermalization as defined by equation ([7|; 
T-(*'')Tg*^ is the thermalization time in units of the initial effective temperature; T^fJ/T^tr is the ratio of the effective 
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VI B 



temperature at thermalization to the initial effective temperature; s„_g„ is dimensionless initial entropy density (77) 

defined by apparent horizon, as described in section 
fluid hydrodynamics. 
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Appendix B: Full Einstein equations 

The following are Einstein equations for all ADM functions. Due to the size of expressions for L, M, P, they were 
exported from a notebook and follow standard Mathematica notation for partial derivatives. 
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Equation for the evolution of lapse function a(t,u) was obtained from the algebraic expression specifying it by 
differentiation with respect to time. 



[1] U. W. Heinz, "Thermalization at RHIC," AIP Conf. 
Proc. 739, 163 (2005) [arXiv:nucl-th/0407067] . 

[2] P. M. Chesler and L. G. Yaffe, "Horizon formation 
and far-from-equilibrium isotropization in supersymmet- 
ric Yang-Mills plasma," Phys. Rev. Lett. 102, 211601 



(2009) [arXiv:0812.2053 [hep-th]]. 
[3] M. P. Heller, R. A. Janik, P. Witaszczyk, "The char- 
acteristics of thermalization of boost-invariant plasma 
from holography," Phys. Rev. Lett. 108, 201602 (2012) 
[arXiv: 1103.3452 [hep-th]]. 



23 



[4] Michael Strickland's talk "Is early isotropization in heavy 
ion c ollisions a necessary condition to describe HIC 
data?" at EMMI Rapid Reaction Task Force on Ther- 



[22; 



[5 

[6: 

[7; 

[8 
[9 

[10: 

[11 

[12 
[13 
[14 
[15 
[16 
[17 
[18 
[19 
[20; 

[21 



malization in Nonabelian Plasmas, University of Heidel- [23 
berg, Germany, 2011. For the conference summary, see 

m- 

J. Berges, J. -P. Blaizot and F. Gelis, "EMMI Rapid 
Reaction Task Force on 'Thermalization in Non-abelian 
Plasmas'," arXiv: 1203.2042 [hep-ph]. [24 

S. Bhattacharyya, R. Loganayagam, I. Mandal, S. Min- 
walla and A. Sharma, "Conformal Nonlinear Fluid Dy- 
namics from Gravity in Arbitrary Dimensions," JHEP [25 
0812, 116 (2008) [arXiv:0809.4272 [hep-th]]. 
G. T. Horowitz and V. E. Hubeny, "Quasinormal Modes [26 
of AdS Black Holes and the Approach to Thermal Equi- 
hbrium," Phys. Rev. D 62, 024027 (2000) [arXiv:hep- 
th/9909056v2]. 

R. Janik and R. Peschanski, "Gauge/gravity duality and [27 
thermalization of a boost-invariant perfect fluid," Phys. 
Rev. D 74, 046007 (2006) [arXiv:hep-th/0606149v3]. 
M. P. Heller, D. Mateos, W. van der Schee and D. Tran- [28 
canelli, "Strong Coupling Isotropization of Non-Abelian 
Plasmas Simplified," Phys. Rev. Lett. 108, 191601 (2012) 
[arXiv:1202.0981 [hep-th]]. 

P. M. Chesler and L. G. Yaffe, "Boost invariant flow, [29 
black hole formation, and far-from-equilibrium dynamics 
in N = 4 supersymmetric Yang-Mills theory," Phys. Rev. 
D 82, 026006 (2010) [arXiv:0906.4426 [hep-th]]. [30 

P. M. Chesler and L. G. Yaffe, "Holography and 
colliding gravitational shock waves in asymptotically 
AdSs spacetime," Phys. Rev. Lett. 106, 021601 (2011) [31 

[arXiv:1011.3562 [hep-th]]. 

P. Bizon and A. Rostworowski, "On weakly turbulent in- 
stability of anti-de Sitter space," Phys. Rev. Lett 107, 
031102 (2011) [arXiv:1104.3702 [gr-qc]]. [32 

0. J. C. Dias, G. T. Horowitz and J. E. Santos, "Grav- 
itational Turbulent Instability of Anti-de Sitter Space," 
[arXiv:1109.1825vl [hep-th]]. [33 
D. Garfinkle, L. A. P. Zayas, "Rapid Thermaliza- 
tion in Field Theory from Gravitational Collapse," 
arXiv:1106.2339v3. [34 
S. Bhattacharyya, V. E. Hubeny, S. Minwalla, M. Ranga- 
mani, "Nonlinear Fluid Dynamics from Gravity," JHEP 
0802, 045 (2008). [arXiv:0712.2456 [hep-th]]. [35 
J. D. Bjorken, "Highly Relativistic Nucleus-Nucleus Col- 
lisions: The Central Rapidity Region," Phys. Rev. D 27, [36 
140 (1983). 

R. A. Janik and R. B. Peschanski, "Asymptotic perfect 
fluid dynamics as a consequence of AdS/CFT," Phys. [37 
Rev. D 73, 045013 (2006) [arXiv:hep-th/0512162]. 
R. A. Janik, "Viscous plasma evolution from gravity 
using AdS/CFT," Phys. Rev. Lett. 98, 022302 (2007) [38 

[arXiv:hep-th/0610144]. 

M. P. Heller and R. A. Janik, "Viscous hydrodynam- 
ics relaxation time from AdS/CFT," Phys. Rev. D 76, [39 
025027 (2007) [arXiv;hep-th/0703243]. 

1. Booth, M. P. Heller and M. Spalinski, "Black brane [40 
entropy and hydrodynamics: the boost-invariant case," 
Phys. Rev. D 80, 126013 (2009) [arXiv:0910.0748 [hep- 
th]]. [41 
M. P. Heller, P. Surowka, R. Loganayagam, M. Spalinski, 

S. E. Vazquez, "Consistent Holographic Description of 
Boost-Invariant Plasma," Phys. Rev. Lett. 102, 041601 [42 
(2009). [arXiv:0805.3774 [hep-th]]. 



S. Kinoshita, S. Mukohyama, S. Nakamura and 
K. y. Oda, "A Holographic Dual of Bjorken Flow," Prog. 
Theor. Phys. 121, 121 (2009) [arXiv:0807.3797 [hep-th]]. 
S. Kinoshita, S. Mukohyama, S. Nakamura and 
K. y. Oda, "Consistent Anti-de Sitter-Space/Conformal- 
Field-Theory Dual for a Time-Dependent Finite Tem- 
perature System," Phys. Rev. Lett. 102, 031601 (2009) 
[arXiv:0901.4834 [hep-th]]. 

G. Beuf, M. P. Heller, R. A. Janik and R. Peschanski, 
"Boost-invariant early time dynamics from AdS/CFT," 
JHEP 0910, 043 (2009) [arXiv:0906.4423 [hep-th]]. 
C. W. Misner, K. S. Thorne, J. A. Wheeler, "Gravita- 
tion," San Francisco 1973, 1279p. 

H. Bantilan, F. Pretorius and S. S. Gubser, "Simula- 
tion of Asymptotically AdS5 Spacetimes with a Gen- 
eralized Harmonic Evolution Scheme," [arXiv:1201.2132 
[hep-th]]. 

H. Yoshino and M. Shibata, "Higher-dimensional numer- 
ical relativity: Formulation and code tests," Phys. Rev. 
D 80, 084025 (2009) [arXiv:0907.2760v2 [gr-qc]]. 
M. Lublinsky and E. Shuryak, "How much entropy is pro- 
duced in strongly coupled Quark-Gluon Plasma (sQGP) 
by dissipative effects?," Phys. Rev. C 76, 021901 (2007) 
[arXiv:0704.1647 [hep-ph]]. 

M. Lublinsky and E. Shuryak, "Improved Hydrodynam- 
ics from the AdS/CFT," Phys. Rev. D 80, 065026 (2009) 
[arXiv:0905.4069 [hep-ph]]. 

A. Buchel, "Shear viscosity of boost invariant plasma 
at flnite couphng," Nucl. Phys. B 802, 281 (2008) 
[arXiv:0801.4421 [hep-th]]. 

S. de Haro, S. N. Solodukhin, K. Skenderis, "Holographic 
reconstruction of space-time and renormalization in the 
AdS / GET correspondence," Commun. Math. Phys. 
217, 595-622 (2001). [hep-th/0002230]. 
E. Poisson, "A Relativist's Toolkit: The Mathematics 
of Black-Hole Mechanics," Cambridge University Press 
(2004). 

M. Alcubierre, "Hyperbolic slicings of space-time: Singu- 
larity avoidance and gauge shocks," Class. Quant. Grav. 
20, 607 (2003) [gr-qc/02 10050]. 

T. W. Baumgarte and S. L. Shapiro, "Numerical relativ- 
ity and compact binaries," Phys. Rept. 376, 41 (2003) 
[gr-qc/0211028]. 

E. Gourgoulhon, "3+1 formalism and bases of numerical 
relativity," gr-qc/0703035 [GR-QC]. 

P. Figueras, V. E. Hubeny, M. Rangamani and S. F. Ross, 
"Dynamical black holes and expanding plasmas," JHEP 
0904, 137 (2009) [arXiv:0902.4696 [hep-th]]. 
S. Bhattacharyya et al., "Local Fluid Dynamical Entropy 
from Gravity,'' JHEP 0806, 055 (2008) [arXiv:0803.2526 
[hep-th]]. 

A. Ashtekar and B. Krishnan, "Isolated and dynamical 
horizons and their applications," Living Rev. Rel. 7, 10 

(2004) [gr-qc/0407042]. 

I. Booth, "Black hole boundaries," Can. J. Phys. 83, 1073 

(2005) [arXiv:gr-qc/0508107]. 

I. Booth, M. P. Heller, G. Plewa and M. Spalinski, "On 
the apparent horizon in fluid-gravity duality," Phys. Rev. 
D 83, 106005 (2011) [arXiv:1102.2885 [hep-th]]. 
V. E. Hubeny, M. Rangamani and T. Takayanagi, "A 
Covariant holographic entanglement entropy proposal," 
JHEP 0707, 062 (2007) [arXiv:0705.0016 [hep-th]]. 
V. Balasubramanian and P. Kraus, "A Stress tensor for 
Anti-de Sitter gravity," Commun. Math. Phys. 208, 413 



24 



(1999) [hep-th/9902121]. tion," arXiv:1103.2683 [hep-th]. 

[43] J. Abajo-Arrastia, J. Aparicio and E. Lopez, "Holo- [46] P. Romatschke and R. Venugopalan, "Collective non- 
graphic Evolution of Entanglement Entropy," JHEP Abelian instabilities in a melting color glass condensate," 
1011, 149 (2010) [arXiv: 1006.4090 [hep-th]]. Phys. Rev. Lett. 96, 062302 (2006) [hep-ph/0510121]. 

[44] V. Balasubramanian et al., "Thormalization of Strongly [47] R. Snellings, "Elliptic Flow: A Brief Review," New J. 
Coupled Field Theories," Phys. Rev. Lett. 106, 191601 Phys. 13, 055008 (2011) [arXiv: 1102.3010 [nucl-ex]]. 

(2011) [arXiv: 1012.4753 [hep-th]]. [48] B. MuUcr and A. Schafer, "Transverse energy density 

[45] V. Balasubramanian et al., "Holographic Thermaliza- fluctuations in the Color Glass Condensate Model," 

arXiv:1111.3347 [hep-ph]. 



25 



